Implement support for stoichiometries with the same name.

doc/AntimonyTutorial.htm

@@ -755,14 +755,47 @@ <h3 id="substance-only-species">Substance-only species</h3>
 <p>Now, whenever ‘S1’ is used in the model, it is a reference to the
 species amount, and not its concentration. Defining an initial amount is
 also changed:</p>
-<p>S1 = 2.5;</p>
+<pre><code>S1 = 2.5;</code></pre>
 <p>This will set the initial amount to 2.5, not the initial
 concentration. If you wish to set the initial concentration instead, use
 the compartment:</p>
-<p>S1 = 3.1*C</p>
+<pre><code>S1 = 3.1*C</code></pre>
 <p>Because a concentration times the compartment volume yields an
 amount, in this formulation, ‘3.1’ is set as the initial
 concentration.</p>
+<h3 id="named-stoichiometries">Named stoichiometries</h3>
+<p>A stoichiometry in a reaction may be given an ID instead of a number,
+and that ID may be set later:</p>
+<pre><code>J0: n A -&gt; B; k1*A^n
+n = 3</code></pre>
+<p>The id of the stoichiometry may now be changed by other model
+constructs: events, rate rules, and assignment rules may all use the
+value as a target:</p>
+<pre><code>J0: n A -&gt; m B; k1*A^n
+n := time/3
+m = 1
+at A &lt; 3: m = 2</code></pre>
+<p>This also gives the stoichiometry an ID that can be given a value
+directly (or be tracked) by some simulators (such as roadrunner).</p>
+<p>If you want to use the same ID for multiple stoichiometries, this can
+be done straightforwardly:</p>
+<pre><code>J0: n A -&gt; n B; k1*A^n
+n = 3</code></pre>
+<p>However! When translated to SBML, every stoichiometry must have a
+unique ID. Therefore, an assignment rule will be created to set the
+value of B’s stoichiometry to ‘n’. Effectively, the model will
+become:</p>
+<pre><code>J0: n A -&gt; J0_B_stoich B; k1*A^n
+n = 3
+J0_B_stoich := n</code></pre>
+<p>This is mathematically identical, but some simulators may balk at an
+assignment rule to a stoichiometry, as this is a feature of SBML that
+not everyone supports. If this happens, just name all your
+stoichiometries uniquely:</p>
+<pre><code>J0: n A -&gt; m B; k1*A^n
+n = 3
+m = 3</code></pre>
+<p>and remember to change them both at the same time.</p>
 <h3 id="modules">Modules</h3>
 <p>Antimony input files may define several different models, and may use
 previously-defined models as parts of newly-defined models. Each

doc/AntimonyTutorial.md

@@ -850,16 +850,51 @@ Now, whenever ‘S1’ is used in the model, it is a reference to the
 species amount, and not its concentration. Defining an initial amount is
 also changed:
 
-S1 = 2.5;
+    S1 = 2.5;
 
 This will set the initial amount to 2.5, not the initial concentration.
 If you wish to set the initial concentration instead, use the compartment:
 
-S1 = 3.1*C
+    S1 = 3.1*C
 
 Because a concentration times the compartment volume yields an amount, in
 this formulation, '3.1' is set as the initial concentration.
 
+### Named stoichiometries
+
+A stoichiometry in a reaction may be given an ID instead of a number, and that ID may be set later:
+
+    J0: n A -> B; k1*A^n
+    n = 3
+
+The id of the stoichiometry may now be changed by other model constructs:  events, rate rules, and assignment rules may all use the value as a target:
+
+    J0: n A -> m B; k1*A^n
+    n := time/3
+    m = 1
+    at A < 3: m = 2
+
+This also gives the stoichiometry an ID that can be given a value directly (or be tracked) by some simulators (such as roadrunner).
+
+If you want to use the same ID for multiple stoichiometries, this can be done straightforwardly:
+
+    J0: n A -> n B; k1*A^n
+    n = 3
+
+However!  When translated to SBML, every stoichiometry must have a unique ID.  Therefore, an assignment rule will be created to set the value of B's stoichiometry to 'n'.  Effectively, the model will become:
+
+    J0: n A -> J0_B_stoich B; k1*A^n
+    n = 3
+    J0_B_stoich := n
+
+This is mathematically identical, but some simulators may balk at an assignment rule to a stoichiometry, as this is a feature of SBML that not everyone supports.  If this happens, just name all your stoichiometries uniquely:
+
+    J0: n A -> m B; k1*A^n
+    n = 3
+    m = 3
+
+and remember to change them both at the same time.
+
 ### Modules
 
 Antimony input files may define several different models, and may use

src/module-sbml.cpp

@@ -1502,19 +1502,52 @@ void Module::LoadSBML(Model* sbml)
           stoichvar->SetAssignmentRule(&formula);
         }
         else if (specref->isSetIdAttribute()) {
-          stoichvar = AddOrFindVariable(&(specref->getIdAttribute()));
-          bool setret = stoichvar->SetType(varStoichiometry); //Since the SBML file is valid.
-          assert(!setret);
-          if (specref->isSetStoichiometry() && !stoichvar->HasFormula()) {
-            Formula formula;
-            formula.AddNum(specref->getStoichiometry());
-            if (specref->getLevel() == 1 && specref->getDenominator() != 1) {
-              formula.AddMathThing('/');
-              formula.AddNum(specref->getDenominator());
+          string stoichid = specref->getIdAttribute();
+          //Restrictive check: only collapse this speciesReference back to a
+          //shared Antimony symbol (rather than giving it its own named
+          //stoichiometry) if it exactly matches what our writer generates
+          //for a reused stoichiometry: SBO:0000481 ("stoichiometric
+          //coefficient"), an id of the form "<reaction>_<species>_stoich[N]",
+          //and an assignment rule that's nothing but a bare reference to
+          //another stoichiometry.
+          Variable* collapseTarget = NULL;
+          if (specref->getSBOTerm() == 481) {
+            string expectedPrefix = reactionName + "_" + specref->getSpecies() + "_stoich";
+            if (stoichid.rfind(expectedPrefix, 0) == 0) {
+              string suffix = stoichid.substr(expectedPrefix.size());
+              bool nameMatches = suffix.empty() ||
+                  suffix.find_first_not_of("0123456789") == string::npos;
+              if (nameMatches) {
+                const Rule* aliasrule = sbml->getRule(stoichid);
+                if (aliasrule != NULL && aliasrule->isAssignment() && aliasrule->isSetMath() &&
+                    aliasrule->getMath()->getType() == AST_NAME) {
+                  string targetname = aliasrule->getMath()->getName();
+                  Variable* target = AddOrFindVariable(&targetname);
+                  if (target->GetType() == varStoichiometry) {
+                    collapseTarget = target;
+                  }
+                }
+              }
+            }
+          }
+          if (collapseTarget != NULL) {
+            stoichvar = collapseTarget;
+          }
+          else {
+            stoichvar = AddOrFindVariable(&stoichid);
+            bool setret = stoichvar->SetType(varStoichiometry); //Since the SBML file is valid.
+            assert(!setret);
+            if (specref->isSetStoichiometry() && !stoichvar->HasFormula()) {
+              Formula formula;
+              formula.AddNum(specref->getStoichiometry());
+              if (specref->getLevel() == 1 && specref->getDenominator() != 1) {
+                formula.AddMathThing('/');
+                formula.AddNum(specref->getDenominator());
+              }
+              stoichvar->SetFormula(&formula);
             }
-            stoichvar->SetFormula(&formula);
+            TranslateRulesAndAssignmentsTo(specref, stoichvar);
           }
-          TranslateRulesAndAssignmentsTo(specref, stoichvar);
         }
         else {
           if (specref->isSetStoichiometry()) {
@@ -2538,6 +2571,9 @@ void Module::CreateSBMLModel(bool comp)
 
   //Reactions
   size_t numrxns = GetNumVariablesOfType(allReactions, comp);
+  //Tracks which named stoichiometry variables have already claimed their
+  //own name as an SBML id in this document; reused ones get an alias id.
+  set<const Variable*> usedNamedStoichIds;
   for (size_t rxn=0; rxn < numrxns; rxn++) {
     const Variable* rxnvar = GetNthVariableOfType(allReactions, rxn, comp);
     const AntimonyReaction* reaction = rxnvar->GetReaction();
@@ -2602,7 +2638,9 @@ void Module::CreateSBMLModel(bool comp)
                 sr->setConstant(true);
                 sr->setStoichiometry(nthstoich);
             }
-            else {
+            else if (usedNamedStoichIds.insert(namedstoich).second) {
+                //First time this stoichiometry variable's name is being used
+                //as an SBML id in this document: keep it as-is.
                 sr->setIdAttribute(namedstoich->GetNameDelimitedBy(cc));
                 if (namedstoich->GetFormulaType() != formulaINITIAL) {
                     sr->setConstant(false);
@@ -2615,6 +2653,35 @@ void Module::CreateSBMLModel(bool comp)
                 }
                 SetAssignmentFor(sbmlmod, namedstoich, syncmap, comp, referencedVars);
             }
+            else {
+                //This stoichiometry variable's name has already been claimed
+                //by an earlier SpeciesReference in this document (Antimony
+                //allows the same symbol to be reused as a stoichiometry more
+                //than once, but SBML ids must be unique).  Give this one its
+                //own id, driven by an assignment rule back to the original,
+                //and mark it as an alias with SBO:0000481 ("stoichiometric
+                //coefficient") plus a "<reaction>_<species>_stoich[N]" name,
+                //so that a restrictive SBML importer can recognize and
+                //collapse it back to the shared symbol.
+                string aliasbase = rxnvar->GetNameDelimitedBy("_") + "_" + nthr->GetNameDelimitedBy("_") + "_stoich";
+                string aliasname = aliasbase;
+                int aliassuffix = 1;
+                while (sbmlmod->getElementBySId(aliasname) != NULL) {
+                    aliassuffix++;
+                    stringstream aliasnum;
+                    aliasnum << aliassuffix;
+                    aliasname = aliasbase + aliasnum.str();
+                }
+                Variable aliasvar(aliasname, this);
+                aliasvar.SetType(varStoichiometry);
+                Formula aliasformula;
+                aliasformula.AddVariable(namedstoich);
+                aliasvar.SetAssignmentRule(&aliasformula);
+                sr->setIdAttribute(aliasname);
+                sr->setSBOTerm(481);
+                sr->setConstant(false);
+                SetAssignmentFor(sbmlmod, &aliasvar, syncmap, comp, referencedVars);
+            }
         }
     }
     //Find 'modifiers' and add them.

src/test/TestAntimonyNamedStoich.cpp

@@ -58,3 +58,13 @@ TEST(AntimonyNamedStoich, test_NamedStoichRate_file)
 {
     compareFileTranslation("namedstoich_rate");
 }
+
+TEST(AntimonyNamedStoich, test_stoichiometries_with_same_name)
+{
+    compareFileTranslationWithDifferences("stoichiometries_with_same_name");
+}
+
+TEST(AntimonyNamedStoich, test_stoichiometries_with_same_name_in_multiple_reactions)
+{
+    compareFileTranslationWithDifferences("stoichiometries_with_same_name_in_multiple_reactions");
+}

src/test/TestAntimonyUtil.cpp

@@ -12,6 +12,8 @@
 #include "stringx.h"
 
 #include <string>
+#include <fstream>
+#include <sstream>
 #include "gtest/gtest.h"
 
 using namespace std;
@@ -49,6 +51,39 @@ void compareFileTranslation(const string& base)
   freeAll();
 }
 
+void compareFileTranslationWithDifferences(const string& base)
+{
+  clearPreviousLoads();
+  // load document
+  string dir(TestDataDirectory);
+  string filename = dir + base + ".txt";
+  long ret = loadAntimonyFile(filename.c_str());
+  EXPECT_TRUE(ret != -1);
+  char* atosbml = getCompSBMLString(NULL);
+  EXPECT_TRUE(atosbml != NULL);
+
+  string sbmlfile = dir + base + ".xml";
+  SBMLDocument* doc = readSBMLFromFile(sbmlfile.c_str());
+  string matching = writeSBMLToStdString(doc);
+  EXPECT_STREQ(atosbml, matching.c_str());
+
+  //Now check the roundtripped version, comparing directly against the
+  //reference file's literal contents instead of reparsing it through
+  //Antimony (which would reorder declarations to match direct-parse order).
+  ret = loadSBMLString(matching.c_str());
+  char* roundtrip = getAntimonyString(NULL);
+  EXPECT_TRUE(roundtrip != NULL);
+
+  string rtfilename = dir + base + "_rt.txt";
+  std::ifstream t(rtfilename.c_str());
+  std::stringstream rtref;
+  rtref << t.rdbuf();
+  EXPECT_STREQ(NormalizeLineEndings(rtref.str()).c_str(), NormalizeLineEndings(string(roundtrip)).c_str());
+
+  delete doc;
+  freeAll();
+}
+
 void compareFileTranslationWithRenaming(const string& base)
 {
   clearPreviousLoads();

src/test/TestAntimonyUtil.h

@@ -25,6 +25,18 @@ void compareFileTranslation(const std::string& base);
 // <base>.txt, then converts back to SBML and confirms it round-trips.
 void compareFileTranslationWithRenaming(const std::string& base);
 
+// Like compareFileTranslation, but for the final round-trip check, compares
+// the SBML-loaded antimony text directly against the literal contents of
+// <base>_rt.txt instead of reparsing that file through Antimony and
+// reprinting it. Use this when the SBML path and a direct top-to-bottom
+// Antimony parse are expected to produce the same content but in a
+// different variable declaration order (e.g. because SBML visits
+// Parameters before Reactions, while Antimony declares symbols in the
+// order they're first mentioned in a reaction), so reparsing the reference
+// file would just reintroduce the direct-parse order instead of confirming
+// a match.
+void compareFileTranslationWithDifferences(const std::string& base);
+
 // Loads the given antimony string, compares its SBML translation against
 // the SBML in the given file (relative to TestDataDirectory).
 void compareStringTranslation(const std::string& antimony, const std::string& sbml);

src/test/test-data/stoichiometries_with_same_name.txt

@@ -0,0 +1 @@
+J0: n A -> n B; k1*A^n

src/test/test-data/stoichiometries_with_same_name.xml

@@ -0,0 +1,46 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!-- Created by libAntimony version v3.2.0 with libSBML version 5.21.1. -->
+<sbml xmlns="http://www.sbml.org/sbml/level3/version2/core" xmlns:comp="http://www.sbml.org/sbml/level3/version1/comp/version1" level="3" version="2" comp:required="true">
+  <model metaid="__main" id="__main">
+    <listOfCompartments>
+      <compartment sboTerm="SBO:0000410" id="default_compartment" spatialDimensions="3" size="1" constant="true"/>
+    </listOfCompartments>
+    <listOfSpecies>
+      <species id="A" compartment="default_compartment" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
+      <species id="B" compartment="default_compartment" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
+    </listOfSpecies>
+    <listOfParameters>
+      <parameter id="k1" constant="true"/>
+    </listOfParameters>
+    <listOfRules>
+      <assignmentRule variable="J0_B_stoich">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <ci> n </ci>
+        </math>
+      </assignmentRule>
+    </listOfRules>
+    <listOfReactions>
+      <reaction id="J0" reversible="true">
+        <listOfReactants>
+          <speciesReference id="n" species="A" constant="true"/>
+        </listOfReactants>
+        <listOfProducts>
+          <speciesReference sboTerm="SBO:0000481" id="J0_B_stoich" species="B" constant="false"/>
+        </listOfProducts>
+        <kineticLaw>
+          <math xmlns="http://www.w3.org/1998/Math/MathML">
+            <apply>
+              <times/>
+              <ci> k1 </ci>
+              <apply>
+                <power/>
+                <ci> A </ci>
+                <ci> n </ci>
+              </apply>
+            </apply>
+          </math>
+        </kineticLaw>
+      </reaction>
+    </listOfReactions>
+  </model>
+</sbml>

src/test/test-data/stoichiometries_with_same_name_in_multiple_reactions.txt

@@ -0,0 +1,3 @@
+J0: n A -> n B; k1*A^n
+J1: n C -> n B; k1*C^n
+J2: m A + n B -> m C + n D; k2