fix: preserve deterministic behavior with connected updates

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; r<firedReaction->getNumConnectedRxns(); r++) {
-		rxn = firedReaction->getconnectedRxn(r);
-		for (int pos=0; pos<rxn->getNumOfReactants(); 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; r<reactions.size(); r++) {
+		rxn = reactions.at(r);
+		if (!this->system->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
 // }
 
 
-

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; i<allReactantTemplates.size(); i++) {
+		if (rxn->isTemplateCompatible(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; i<allProductTemplates.size(); i++) {
+			if (rxn->isTemplateCompatible(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<Molecule*> directProductSet;
+	for (unsigned int msIndex=0; msIndex<n_mappingsets; msIndex++) {
+		MappingSet *ms = mappingSet[msIndex];
+		if (ms==0) continue;
+		for (unsigned int mapIndex=0; mapIndex<ms->getNumOfMappings(); 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

src/NFcore/templateMolecule.cpp

@@ -2120,4 +2120,4 @@ bool TemplateMolecule::isMoleculeTypeAndComponentPresent(MoleculeType * mt, int
 	}
 
 	return false;
-}
+}

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; r<n_reactants; r++) {
 		for (unsigned int i=0; i<transformations[r].size(); i++) {
 			transfn = transformations[r].at(i);
 			t1 = transfn->getTemplateMolecule();
 			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;
 		}

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: