diff --git a/src/NFcore/NFcore.hh b/src/NFcore/NFcore.hh index 95b4dfe4..7aad6297 100644 --- a/src/NFcore/NFcore.hh +++ b/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 complexMembers; list ::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: }; diff --git a/src/NFcore/complex.cpp b/src/NFcore/complex.cpp index f38b66ab..b9ca09d6 100644 --- a/src/NFcore/complex.cpp +++ b/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!!! "< 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: " <getComplexID()<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(); + } } } diff --git a/src/NFcore/reactionClass.cpp b/src/NFcore/reactionClass.cpp index 6e11c9f8..b669d4fb 100755 --- a/src/NFcore/reactionClass.cpp +++ b/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 updatedComplexIds; diff --git a/src/NFcore/system.cpp b/src/NFcore/system.cpp index 0393d936..0dbb7439 100644 --- a/src/NFcore/system.cpp +++ b/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; kstraightAdd(); - } + complex->ensureSpeciesObsCache(nSpeciesObs); + + if (complex->isSpeciesObsDirty()) { + for (int i=0; iisObservable(complex); + complex->getSpeciesObsCache()[i] = match; + } + complex->clearSpeciesObsDirty(); + } + + for (int i=0; igetSpeciesObsCache()[i]; + for (int k=0; kstraightAdd(); + } } } - /* - 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; kstraightAdd(); - } - } - } - */ } diff --git a/src/NFinput/NFinput.cpp b/src/NFinput/NFinput.cpp index 943fdbf3..ea89f165 100644 --- a/src/NFinput/NFinput.cpp +++ b/src/NFinput/NFinput.cpp @@ -1746,14 +1746,83 @@ bool NFinput::initReactionRulePermutation( } } - if (pRxnRule->FirstChildElement("ListOfExcludeProducts") || - pRxnRule->FirstChildElement("ListOfIncludeProducts")) + // Parse ListOfExcludeProducts + // Format: contains 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 dummyComps, dummySymMap; + map 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: contains 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 dummyComps, dummySymMap; + map 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 diff --git a/src/NFreactions/transformations/transformationSet.cpp b/src/NFreactions/transformations/transformationSet.cpp index 22df5684..734a8b89 100644 --- a/src/NFreactions/transformations/transformationSet.cpp +++ b/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& 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& parsedTemplates) { + ProductFilter pf; + pf.productIndex = productIndex; + pf.pattern = pattern; + pf.isExclude = false; + pf.parsedTemplates = parsedTemplates; + productFilters.push_back(pf); +} + +bool TransformationSet::checkProductFilters(const list &products) const { + if (productFilters.empty()) return true; + + // Collect unique complexes from the product molecule list + unordered_set 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; +} diff --git a/src/NFreactions/transformations/transformationSet.hh b/src/NFreactions/transformations/transformationSet.hh index b690208d..524615a0 100644 --- a/src/NFreactions/transformations/transformationSet.hh +++ b/src/NFreactions/transformations/transformationSet.hh @@ -50,6 +50,13 @@ namespace NFcore map parsedTemplates; }; + struct ProductFilter { + int productIndex; + TemplateMolecule *pattern; + bool isExclude; + map parsedTemplates; + }; + public: /*! @@ -315,6 +322,10 @@ namespace NFcore void addIncludeReactant(int reactantIndex, TemplateMolecule *pattern, const map& parsedTemplates); bool checkReactantFilters(int reactantIndex, Molecule *mol) const; + void addExcludeProduct(int productIndex, TemplateMolecule *pattern, const map& parsedTemplates); + void addIncludeProduct(int productIndex, TemplateMolecule *pattern, const map& parsedTemplates); + bool checkProductFilters(const list &products) const; + protected: bool addBindingTransformImpl(TemplateMolecule *t1, string bSiteName1, TemplateMolecule *t2, string bSiteName2, bool isNewMolecule); @@ -395,6 +406,7 @@ namespace NFcore vector < pair > collision_pairs; vector reactantFilters; + vector productFilters; private: void initCommon();