diff --git a/src/main/java/fern/example/SBMLMathTreeTest.java b/src/main/java/fern/example/SBMLMathTreeTest.java deleted file mode 100644 index 9fa4628b..00000000 --- a/src/main/java/fern/example/SBMLMathTreeTest.java +++ /dev/null @@ -1,32 +0,0 @@ -package fern.example; - -import fern.network.FeatureNotSupportedException; -import fern.network.Network; -import fern.network.sbml.SBMLNetwork; -import fern.network.sbml.SBMLPropensityCalculator; -import fern.tools.NetworkTools; -import org.sbml.jsbml.validator.ModelOverdeterminedException; - -import javax.xml.stream.XMLStreamException; -import java.io.IOException; - -public class SBMLMathTreeTest { - - /** - * Dumb the MathTrees of an SBML network to Stdout - * - * @param args - * @throws FeatureNotSupportedException - */ - public static void main(String[] args) - throws FeatureNotSupportedException, IOException, XMLStreamException, ModelOverdeterminedException { - Network net = new SBMLNetwork(ExamplePath.find("mapk_sbml.xml")); - for (int i = 0; i < net.getNumReactions(); i++) { - System.out.println("Reaction " + NetworkTools.getReactionNameWithAmounts(net, i)); - NetworkTools - .dumpMathTree(((SBMLPropensityCalculator) net.getPropensityCalculator()).getMathTree(i)); - System.out.println(); - } - } - -} diff --git a/src/main/java/fern/network/sbml/MathTree.java b/src/main/java/fern/network/sbml/MathTree.java deleted file mode 100644 index 571abbfc..00000000 --- a/src/main/java/fern/network/sbml/MathTree.java +++ /dev/null @@ -1,88 +0,0 @@ -package fern.network.sbml; - -import fern.network.AmountManager; -import fern.simulation.Simulator; -import java.util.LinkedList; -import java.util.List; -import java.util.Stack; -import org.sbml.jsbml.ASTNode; -import org.simulator.sbml.SBMLinterpreter; -import org.simulator.sbml.astnode.ASTNodeValue; - -/** - * Representation of am evaluation tree. Within a sbml file, MathML branches may occur at different - * positions. These are represented as MathTrees in FERN. - * - * @author Florian Erhard - */ -public class MathTree { - - private ASTNode copiedAST; - private SBMLinterpreter sbmlInterpreter; - public static final String TEMP_VALUE = "SBML_SIMULATION_TEMP_VALUE"; - - /** - * Creates a MathTree {@link ASTNode}. - * - * @param interpreter sbmlInterpreter instance for calculating the nodes - * @param ast ASTNode - */ - public MathTree(SBMLinterpreter interpreter, ASTNode ast) { - sbmlInterpreter = interpreter; - copiedAST = interpreter.copyAST(ast, true, null, null); - } - - /** - * Gets the species present in this tree. - * - * @return indices of the species. - */ - public List getSpecies() { - List re = new LinkedList<>(); - Stack dfs = new Stack<>(); - dfs.add(copiedAST); - while (!dfs.empty()) { - ASTNode node = dfs.pop(); - if ((node.getNumChildren() == 0) && !node.isOperator() && !node.isNumber()) { - Integer index = null; - if (sbmlInterpreter.getModel().getSpecies(node.getName()) != null) { - // Subtracting from the total compartment count as species indices start after compartments - // in the Y array in the interpreter. - index = sbmlInterpreter.getSymbolHash().get(node.getName()) - sbmlInterpreter.getModel() - .getCompartmentCount(); - } - if ((index != null) && !re.contains(index)) { - re.add(index); - } - } else if (node.getNumChildren() != 0) { - for (int i = 0; i < node.getNumChildren(); i++) { - dfs.add(node.getChild(i)); - } - } - } - return re; - } - - /** - * Gets the ASTNode of this MathTree. - * - * @return the copiedAST of this MathTree - */ - public ASTNode getCopiedAST() { - return copiedAST; - } - - /** - * Evaluate the MathTree. - * - * @param amount AmountManager - * @param sim Simulator - * @return value of the expression - */ - public double calculate(AmountManager amount, Simulator sim) { - sbmlInterpreter.updateSpeciesConcentration(amount); - sbmlInterpreter.setCurrentTime(sim.getTime()); - return ((ASTNodeValue) copiedAST.getUserObject(TEMP_VALUE)).compileDouble(sim.getTime(), 0d); - } - -} diff --git a/src/main/java/fern/network/sbml/SBMLEventHandlerObserver.java b/src/main/java/fern/network/sbml/SBMLEventHandlerObserver.java index bc90ff3b..ad3cd7e5 100644 --- a/src/main/java/fern/network/sbml/SBMLEventHandlerObserver.java +++ b/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 variableAssignment; - private Map parameterAssignment; + private SBMLinterpreter interpreter; + private Map variableAssignment; + private Map 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 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; } - -} +} \ No newline at end of file diff --git a/src/main/java/fern/network/sbml/SBMLPropensityCalculator.java b/src/main/java/fern/network/sbml/SBMLPropensityCalculator.java index 267c95ee..620b2297 100644 --- a/src/main/java/fern/network/sbml/SBMLPropensityCalculator.java +++ b/src/main/java/fern/network/sbml/SBMLPropensityCalculator.java @@ -1,42 +1,33 @@ -/* - * 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 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()); @@ -44,52 +35,54 @@ public SBMLPropensityCalculator(SBMLinterpreter interpreter) throws ModelOverdet 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 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 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 getKineticLawSpecies(int reaction) { - return propensities[reaction].getSpecies(); + public List getKineticLawSpecies(int reactionIndex) { + // Instead of doing a manual DFS tree search, we just ask the SBML Reaction for its dependencies + List 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; } - - -} +} \ No newline at end of file diff --git a/src/main/java/fern/tools/NetworkTools.java b/src/main/java/fern/tools/NetworkTools.java index 2bb40d00..5d4ffd04 100644 --- a/src/main/java/fern/tools/NetworkTools.java +++ b/src/main/java/fern/tools/NetworkTools.java @@ -15,7 +15,6 @@ import fern.network.KineticConstantPropensityCalculator; import fern.network.Network; import fern.network.NetworkLoader; -import fern.network.sbml.MathTree; import fern.simulation.Simulator; import fern.simulation.algorithm.GillespieSimple; import org.sbml.jsbml.ASTNode; @@ -355,40 +354,6 @@ public static void dumpNetwork(Network net, Writer writer) throws IOException { writer.flush(); } - /** - * Dumps the MathTree to stdout. - * - * @param tree mathtree - */ - public static void dumpMathTree(MathTree tree) { - try { - dumpMathTree(tree, new PrintWriter(System.out)); - } catch (IOException e) { - e.printStackTrace(); - } - } - - /** - * Dumps the math tree to the given {@link Writer}. - * - * @param tree mathtree - */ - public static void dumpMathTree(MathTree tree, Writer writer) throws IOException { - dumpMathTreeNode(tree.getCopiedAST(), writer, new StringBuilder()); - writer.flush(); - } - - private static void dumpMathTreeNode(ASTNode astNode, Writer writer, StringBuilder indend) - throws IOException { - writer.write(indend.toString()); - writer.write(astNode.toString()); - writer.write("\n"); - indend.append(" "); - for (ASTNode child : astNode.getChildren()) { - dumpMathTreeNode(child, writer, indend); - } - indend.deleteCharAt(0); - } /** * Copies the actual amount of the species in the network to its initial amount.