Refactor: Integrate FERN with native SBSCL SBMLinterpreter (Resolves #55)

src/main/java/fern/network/sbml/SBMLEventHandlerObserver.java

@@ -5,52 +5,48 @@
 import fern.simulation.observer.TriggerObserver;
 import java.util.HashMap;
 import java.util.Map;
+import org.sbml.jsbml.ASTNode;
 import org.sbml.jsbml.Event;
 import org.sbml.jsbml.validator.ModelOverdeterminedException;
 import org.simulator.sbml.SBMLinterpreter;
+import org.simulator.sbml.astnode.ASTNodeValue;
 
 /**
  * Observer which handles an event of a sbml model.
+ * Refactored to use native ASTNode evaluation instead of MathTree.
  *
  * @author Florian Erhard
  */
 public class SBMLEventHandlerObserver extends TriggerObserver {
 
   private String name;
-  private MathTree trigger;
-  private MathTree delay;
+  private ASTNode triggerAST;
+  private ASTNode delayAST;
   private SBMLNetwork net;
-  private Map<String, MathTree> variableAssignment;
-  private Map<String, MathTree> parameterAssignment;
+  private SBMLinterpreter interpreter;
+  private Map<String, ASTNode> variableAssignment;
+  private Map<String, ASTNode> parameterAssignment;
   private boolean lastStepTriggered;
 
-  /**
-   * Creates the observer.
-   *
-   * @param sim         the simulator
-   * @param net         the sbml network
-   * @param interpreter the sbmlInterpreter instance to calculate the node values
-   * @param event       the event object of the sbml model
-   */
   public SBMLEventHandlerObserver(Simulator sim, SBMLNetwork net, SBMLinterpreter interpreter,
       Event event) throws ModelOverdeterminedException {
     super(sim);
-
     this.net = net;
+    this.interpreter = interpreter;
     parse(event, interpreter);
   }
 
   private void parse(Event event, SBMLinterpreter interpreter) {
     this.name = event.getId();
-    this.trigger = new MathTree(interpreter, event.getTrigger().getMath());
-    this.delay =
-        event.getDelay() == null ? null : new MathTree(interpreter, event.getDelay().getMath());
+    this.triggerAST = interpreter.copyAST(event.getTrigger().getMath(), true, null, null);
+    this.delayAST = event.getDelay() == null ? null : interpreter.copyAST(event.getDelay().getMath(), true, null, null);
+
     variableAssignment = new HashMap<>();
     parameterAssignment = new HashMap<>();
 
     for (int i = 0; i < event.getNumEventAssignments(); i++) {
       String var = event.getEventAssignment(i).getVariable();
-      MathTree tree = new MathTree(interpreter, event.getEventAssignment(i).getMath());
+      ASTNode tree = interpreter.copyAST(event.getEventAssignment(i).getMath(), true, null, null);
       if (interpreter.getModel().containsSpecies(var)) {
         variableAssignment.put(var, tree);
       } else {
@@ -59,23 +55,27 @@ private void parse(Event event, SBMLinterpreter interpreter) {
     }
   }
 
+  private double evaluate(ASTNode ast, Simulator sim) {
+    interpreter.updateSpeciesConcentration(net.getAmountManager());
+    interpreter.setCurrentTime(sim.getTime());
+    return ((ASTNodeValue) ast.getUserObject("SBML_SIMULATION_TEMP_VALUE")).compileDouble(sim.getTime(), 0d);
+  }
+
   private void executeEvent() {
     for (String var : variableAssignment.keySet()) {
       net.getAmountManager().setAmount(net.getSpeciesByName(var),
-          (long) variableAssignment.get(var).calculate(net.getAmountManager(), getSimulator()));
+          (long) evaluate(variableAssignment.get(var), getSimulator()));
     }
     for (String par : parameterAssignment.keySet()) {
       Map<String, Double> globals = ((SBMLPropensityCalculator) net.getPropensityCalculator())
           .getGlobalParameters();
-      globals
-          .put(par, parameterAssignment.get(par).calculate(net.getAmountManager(), getSimulator()));
+      globals.put(par, evaluate(parameterAssignment.get(par), getSimulator()));
     }
     getSimulator().reinitialize();
   }
 
   @Override
-  public void activateReaction(int mu, double tau, FireType fireType,
-      int times) {
+  public void activateReaction(int mu, double tau, FireType fireType, int times) {
   }
 
   @Override
@@ -93,11 +93,10 @@ public void step() {
 
   @Override
   public boolean trigger() {
-    boolean triggered = trigger.calculate(net.getAmountManager(), getSimulator()) != 0;
+    boolean triggered = evaluate(triggerAST, getSimulator()) != 0;
     if (!lastStepTriggered && triggered) {
       lastStepTriggered = triggered;
-      double delaytime =
-          delay == null ? 0 : delay.calculate(net.getAmountManager(), getSimulator());
+      double delaytime = delayAST == null ? 0 : evaluate(delayAST, getSimulator());
       if (delaytime <= 0) {
         executeEvent();
       } else {
@@ -122,5 +121,4 @@ public void setSimulatorAsync(Simulator sim) {
   public String toString() {
     return name;
   }
-
-}
+}

src/main/java/fern/network/sbml/SBMLPropensityCalculator.java

@@ -1,95 +1,88 @@
-/*
- * Created on 12.03.2007
- *
- * To change the template for this generated file go to
- * Window>Preferences>Java>Code Generation>Code and Comments
- */
 package fern.network.sbml;
 
 import fern.network.AmountManager;
 import fern.network.ComplexDependenciesPropensityCalculator;
 import fern.simulation.Simulator;
+import java.util.ArrayList;
 import java.util.HashMap;
 import java.util.List;
 import java.util.Map;
 import org.sbml.jsbml.Model;
 import org.sbml.jsbml.Reaction;
+import org.sbml.jsbml.SpeciesReference;
+import org.sbml.jsbml.ModifierSpeciesReference;
 import org.sbml.jsbml.validator.ModelOverdeterminedException;
 import org.simulator.sbml.SBMLinterpreter;
 
 /**
- * Propensity calculator which is used for {@link SBMLNetwork}s. The propensities are calculated by
- * using a {@link MathTree} derived by the MathML representation of the kinetic law for each
- * reaction.
- *
- * @author Florian Erhard
+ * Propensity calculator which is used for {@link SBMLNetwork}s. 
+ * Refactored to use native SBSCL SBMLinterpreter for optimized stochastic simulation.
  */
 public class SBMLPropensityCalculator implements ComplexDependenciesPropensityCalculator {
 
-  private MathTree[] propensities;
+  private SBMLinterpreter interpreter;
   private Map<String, Double> globalParameter;
+  private Model model;
 
-  /**
-   * Creates the {@link MathTree}s and parses the parameters.
-   *
-   * @param interpreter instance of the SBMLinterpreter
-   */
   public SBMLPropensityCalculator(SBMLinterpreter interpreter) throws ModelOverdeterminedException {
-
-    Model model = interpreter.getModel();
+    this.interpreter = interpreter;
+    this.model = interpreter.getModel();
+
     globalParameter = new HashMap<>();
     for (int i = 0; i < model.getNumParameters(); i++) {
       globalParameter.put(model.getParameter(i).getId(), model.getParameter(i).getValue());
     }
     for (int i = 0; i < model.getNumCompartments(); i++) {
       globalParameter.put(model.getCompartment(i).getId(), model.getCompartment(i).getSize());
     }
-
-    propensities = new MathTree[model.getNumReactions()];
-
-    for (int i = 0; i < model.getNumReactions(); i++) {
-      Map<String, Double> localParameter = new HashMap<>();
-      Reaction reaction = model.getReaction(i);
-      for (int j = 0; j < reaction.getKineticLaw().getLocalParameterCount(); j++) {
-        localParameter.put(reaction.getKineticLaw().getLocalParameter(j).getId(),
-            reaction.getKineticLaw().getLocalParameter(j).getValue());
-      }
-      propensities[i] = new MathTree(interpreter, reaction.getKineticLaw().getMath());
-    }
+
+    // Notice: We no longer create an array of MathTrees!
   }
 
-  /**
-   * Gets the global parameters.
-   *
-   * @return global parameters
-   */
   public Map<String, Double> getGlobalParameters() {
     return globalParameter;
   }
 
-  public double calculatePropensity(int reaction, AmountManager amount, Simulator sim) {
-    double re = propensities[reaction].calculate(amount, sim);
+  public double calculatePropensity(int reactionIndex, AmountManager amount, Simulator sim) {
+    // 1. Update the SBSCL interpreter with the current stochastic state
+    interpreter.updateSpeciesConcentration(amount);
+    interpreter.setCurrentTime(sim.getTime());
+
+    // 2. Let the highly-optimized SBSCL engine calculate the velocity natively!
+    double re = interpreter.compileReaction(reactionIndex);
+
     if (re < 0) {
       throw new RuntimeException(
-          "The propensity of reaction " + sim.getNet().getReactionName(reaction) + " is negative");
+          "The propensity of reaction " + sim.getNet().getReactionName(reactionIndex) + " is negative");
     }
     return Math.abs(re);
   }
 
-  public List<Integer> getKineticLawSpecies(int reaction) {
-    return propensities[reaction].getSpecies();
+  public List<Integer> getKineticLawSpecies(int reactionIndex) {
+    // Instead of doing a manual DFS tree search, we just ask the SBML Reaction for its dependencies
+    List<Integer> speciesIndices = new ArrayList<>();
+    Reaction reaction = model.getReaction(reactionIndex);
+
+    // Add Reactants
+    for (SpeciesReference reactant : reaction.getListOfReactants()) {
+       int index = getSpeciesIndex(reactant.getSpecies());
+       if (index != -1 && !speciesIndices.contains(index)) speciesIndices.add(index);
+    }
+
+    // Add Modifiers
+    for (ModifierSpeciesReference modifier : reaction.getListOfModifiers()) {
+       int index = getSpeciesIndex(modifier.getSpecies());
+       if (index != -1 && !speciesIndices.contains(index)) speciesIndices.add(index);
+    }
+
+    return speciesIndices;
   }
 
-
-  /**
-   * Gets the internal representation of the sbml kinetic law.
-   *
-   * @param reaction index of the reaction
-   * @return a MathTree representation of the kinetic law
-   */
-  public MathTree getMathTree(int reaction) {
-    return propensities[reaction];
+  // Helper method to resolve species IDs to array indices for FERN
+  private int getSpeciesIndex(String speciesId) {
+     if (interpreter.getSymbolHash().containsKey(speciesId)) {
+        return interpreter.getSymbolHash().get(speciesId) - model.getCompartmentCount();
+     }
+     return -1;
   }
-
-
-}
+}