diff --git a/src/NFcore/moleculeType.cpp b/src/NFcore/moleculeType.cpp index 763de2a4..6a68e30e 100644 --- a/src/NFcore/moleculeType.cpp +++ b/src/NFcore/moleculeType.cpp @@ -572,35 +572,33 @@ void MoleculeType::updateRxnMembership(Molecule * m) void MoleculeType::updateConnectedRxnMembership(Molecule * m, ReactionClass * firedReaction) { - // Replace the iteration over all reactions for the MoleculeType in - // MoleculeType::updateRxnMembership by only the - // connectedReactions for the fired Reaction. This is a much smaller loop - // and skips moleculetypes that are not the TemplateMolecule of the reactant - // in the connected reaction right away. - // Arvind Rasi Subramaniam - // - for (int r=0; rgetNumConnectedRxns(); r++) { - rxn = firedReaction->getconnectedRxn(r); - for (int pos=0; posgetNumOfReactants(); pos++) { - if (rxn->getMoleculeTypeOfReactantTemplate(pos) != this) continue; - double oldA = rxn->get_a(); - double oldAwithTotal = rxn->update_a(); - rxn->tryToAdd(m, pos); - double newA = rxn->update_a(); - this->system->update_A_tot(rxn,oldA,newA); - // Used for debugging to see which reaction rates changed - // upon updating molecule membership - // Arvind Rasi Subramaniam Nov 21, 2018 - if (!this->system->getTrackConnected()) continue; - if (oldAwithTotal != newA) { - this->system->getConnectedRxnFileStream() << - this->system->getGlobalEventCounter() << "\t" << - firedReaction->getName() << "\t" << - m->getMoleculeTypeName() << "\t" << - m->getUniqueID() << "\t" << - rxn->getName() << "\t" << - oldAwithTotal << "\t" << newA << endl; - } + // Preserve the MoleculeType's native reaction order so the connectivity path + // mutates reactant containers in the same sequence as a full membership + // refresh, while still using the precomputed connectivity matrix. + for (unsigned int r=0; rsystem->areReactionsConnected( + firedReaction->getRxnId(), rxn->getRxnId())) { + continue; + } + int pos = reactionPositions.at(r); + double oldA = rxn->get_a(); + double oldAwithTotal = rxn->update_a(); + rxn->tryToAdd(m, pos); + double newA = rxn->update_a(); + this->system->update_A_tot(rxn,oldA,newA); + // Used for debugging to see which reaction rates changed + // upon updating molecule membership + // Arvind Rasi Subramaniam Nov 21, 2018 + if (!this->system->getTrackConnected()) continue; + if (oldAwithTotal != newA) { + this->system->getConnectedRxnFileStream() << + this->system->getGlobalEventCounter() << "\t" << + firedReaction->getName() << "\t" << + m->getMoleculeTypeName() << "\t" << + m->getUniqueID() << "\t" << + rxn->getName() << "\t" << + oldAwithTotal << "\t" << newA << endl; } } } @@ -828,4 +826,3 @@ void MoleculeType::printDetails() const // } - diff --git a/src/NFcore/reactionClass.cpp b/src/NFcore/reactionClass.cpp index 6cc92f0b..b5ed4a0a 100755 --- a/src/NFcore/reactionClass.cpp +++ b/src/NFcore/reactionClass.cpp @@ -273,7 +273,22 @@ void ReactionClass::appendConnectedRxn(ReactionClass * rxn) { bool ReactionClass::isReactionConnected(ReactionClass * rxn) { // First check if any of the operations share MoleculeType and components with // one of the reactant templates of rxn. - return this->transformationSet->checkConnection(rxn); + if (this->transformationSet->checkConnection(rxn)) return true; + + // Full membership refresh revisits every explicit reactant template in the + // fired rule, not only templates that carry direct transformations. + for (unsigned int i=0; iisTemplateCompatible(allReactantTemplates[i])) return true; + } + + // Product templates can also create new compatible mappings, but avoid + // broadening pure-synthesis rules where this over-connects add-only paths. + if (n_reactants > 0) { + for (unsigned int i=0; iisTemplateCompatible(allProductTemplates[i])) return true; + } + } + return false; } ReactionClass::~ReactionClass() @@ -491,6 +506,30 @@ string ReactionClass::fire(double random_A_number, bool track) { // Add newly created molecules to the list of products this->transformationSet->getListOfAddedMolecules(mappingSet,products,traversalLimit); + // Track molecules that were explicitly mapped by this firing. Products added + // through bonded-neighborhood traversal must use the full updater to preserve + // the same membership mutation order as the non-connectivity path. + std::unordered_set directProductSet; + for (unsigned int msIndex=0; msIndexgetNumOfMappings(); mapIndex++) { + Mapping *mapping = ms->get(mapIndex); + if (mapping==0) continue; + Molecule *directMol = mapping->getMolecule(); + if (directMol!=0) directProductSet.insert(directMol); + } + } + bool hasIndirectProducts = false; + for (molIter = products.begin(); molIter != products.end(); molIter++) { + Molecule *mol = *molIter; + if (!mol->isAlive()) continue; + if (directProductSet.find(mol)==directProductSet.end()) { + hasIndirectProducts = true; + break; + } + } + // if complex bookkeeping is on, find all product complexes // (this is useful for updating Species Observables and TypeII functions, so keep the info handy). // NOTE: this is a brute force approach: check complex of each molecule. there may be a more @@ -558,8 +597,13 @@ string ReactionClass::fire(double random_A_number, bool track) { //Update this molcule's reaction membership // NOTE: as a side-effect, DORreactions that depend on molecule-scoped local functions // (typeI relationship) will be updated as long as UTL is set appropriately. - if ( mol->isAlive() ) - mol->updateRxnMembership(this, useConnectivity); + if ( mol->isAlive() ) { + bool useConnectedUpdate = + useConnectivity && + !hasIndirectProducts && + directProductSet.find(mol)!=directProductSet.end(); + mol->updateRxnMembership(this, useConnectedUpdate); + } } // update complex-scoped local functions for typeII dependencies diff --git a/src/NFcore/templateMolecule.cpp b/src/NFcore/templateMolecule.cpp index a958bb2e..ec8ac708 100644 --- a/src/NFcore/templateMolecule.cpp +++ b/src/NFcore/templateMolecule.cpp @@ -2120,4 +2120,4 @@ bool TemplateMolecule::isMoleculeTypeAndComponentPresent(MoleculeType * mt, int } return false; -} \ No newline at end of file +} diff --git a/src/NFreactions/transformations/transformationSet.cpp b/src/NFreactions/transformations/transformationSet.cpp index c1356540..f07891ed 100644 --- a/src/NFreactions/transformations/transformationSet.cpp +++ b/src/NFreactions/transformations/transformationSet.cpp @@ -967,26 +967,20 @@ void TransformationSet::finalize() bool TransformationSet::checkConnection(ReactionClass * rxn) { TemplateMolecule * t1; - MoleculeType * mt1; Transformation * transfn; - int c1; for(unsigned int r=0; rgetTemplateMolecule(); if (!t1) continue; - mt1 = t1->getMoleculeType(); // AS2023 - if this is not a removal, track connections, removal // doesn't give any reaction connections, so skip that if (transfn->getType()!=(int)TransformationFactory::REMOVE) { - c1 = transfn->getComponentIndex(); - // If the moleculetype or component is not present in the other reaction, - // it is not connected - if (!rxn->areMoleculeTypeAndComponentPresent(mt1, c1)) continue; - - // If the TemplateMolecule is 'incompatible' with any of the reactants - // or products, then the reaction is not connected - if (!rxn->isTemplateCompatible(t1)) continue; + bool isCompatible = rxn->isTemplateCompatible(t1); + if (!isCompatible) continue; + // Full membership refresh still removes/re-adds compatible mappings + // even when the changed component is outside the target pattern, + // which can change ReactantList/ReactantTree ordering. // Both checks passed for one op so return true return true; } else { @@ -1003,15 +997,11 @@ bool TransformationSet::checkConnection(ReactionClass * rxn) { if (!t1) continue; t1 = t1->getMappedPartner(); if (!t1) continue; - mt1 = t1->getMoleculeType(); - c1 = transfn->getComponentIndex(); - // If the moleculetype or component is present in the other reaction, - // it is not connected - if (!rxn->areMoleculeTypeAndComponentPresent(mt1, c1)) continue; - - // If the TemplateMolecule is 'incompatible' with any of the reactants - // or products, then the reaction is not connected - if (!rxn->isTemplateCompatible(t1)) continue; + bool isCompatible = rxn->isTemplateCompatible(t1); + if (!isCompatible) continue; + // See note above: compatibility alone is enough to require a + // connected update if the fast path is to preserve full-update + // membership ordering. // Both checks passed for one op so return true return true; } diff --git a/validate/validate.py b/validate/validate.py index 2654f7e2..02686403 100644 --- a/validate/validate.py +++ b/validate/validate.py @@ -6,7 +6,10 @@ import fnmatch import sys import tempfile -import bionetgen +try: + import bionetgen +except ImportError: + bionetgen = None nIterations=30 nfsimPrePath='..' @@ -173,6 +176,9 @@ def _bng_generate(self, outputDirectory, fileNumber): # already generated, no need to rerun BNG return + if bionetgen is None: + self.fail('bionetgen Python package is required to generate XML fixtures') + bngFileName = os.path.join(outputDirectory, 'v{0}.bngl'.format(fileNumber)) bionetgen.run(bngFileName, out=outputDirectory, suppress=True) @@ -204,6 +210,42 @@ def _run_nfsim(self, outputDirectory, fileNumber, runOptions): expect_success=True, ) + def _assert_same_seed_connectivity_parity(self, xmlPath, runOptions, label): + with tempfile.TemporaryDirectory() as tmpdir: + offPath = os.path.join(tmpdir, 'off.gdat') + onPath = os.path.join(tmpdir, 'on.gdat') + + connectOptions = f'{runOptions} -connect'.strip() + self._run_nfsim_xml(xmlPath, offPath, runOptions) + self._run_nfsim_xml(xmlPath, onPath, connectOptions) + + offHeaders, offData = self._load_gdat(offPath) + onHeaders, onData = self._load_gdat(onPath) + + self.assertEqual(offHeaders, onHeaders, f'Connectivity regression changed {label} output columns') + self.assertEqual(offData.shape, onData.shape, f'Connectivity regression changed {label} output shape') + self.assertTrue( + np.array_equal(offData, onData), + f'Connectivity regression changed the same-seed {label} trajectory' + ) + + def test_connectivity_preserves_seeded_tlbr_trajectory(self): + self._assert_same_seed_connectivity_parity( + os.path.join(nfsimPrePath, 'test', 'tlbr', 'tlbr.xml'), + '-sim 1 -oSteps 100 -seed 1', + 'TLBR' + ) + + def test_connectivity_preserves_seeded_local_function_trajectory(self): + # testSuite/t3 exercises local-function membership updates on a much + # smaller model than AN_chemotaxis while still reproducing the + # master-vs-connect divergence fixed by this branch. + self._assert_same_seed_connectivity_parity( + os.path.join(nfsimPrePath, 'test', 'testSuite', 't3.xml'), + '-sim 1 -oSteps 20 -seed 1', + 'testSuite t3' + ) + def _assert_matching_output_schedules(self, xmlName, continuousOptions, chunkedOptions): xmlPath = os.path.join(mfolder, xmlName) with tempfile.TemporaryDirectory(prefix='nfsim_step_to_') as tmpdir: