Fixed some bugs with ODE simulator getting stuck.

This commit is contained in:
leandrohw 2018-03-05 22:53:23 -08:00
parent 05a355587d
commit 8553a3dcef
4 changed files with 73 additions and 49 deletions

View file

@ -494,12 +494,14 @@ public abstract class HierarchicalSimulation extends AbstractSimulator
protected void checkEvents()
{
double time = currentTime.getState().getStateValue();
for(HierarchicalModel model : modules)
{
int index = model.getIndex();
for (EventNode event : model.getEvents())
{
if(event.isTriggeredAtTime(currentTime.getState().getStateValue(), index))
if(event.isTriggeredAtTime(time, index))
{
event.setMaxDisabledTime(index, Double.NEGATIVE_INFINITY);
event.setMinEnabledTime(index, Double.POSITIVE_INFINITY);

View file

@ -76,6 +76,15 @@ public class EventNode extends HierarchicalNode
{
eventState.get(index).maxDisabledTime = maxDisabledTime;
}
public void setInitialTrue(int index, boolean isInitialTrue)
{
eventState.get(index).isInitialTrue = isInitialTrue;
if(!isInitialTrue)
{
eventState.get(index).maxDisabledTime = 0;
}
}
public double getMinEnabledTime(int index)
{
@ -201,6 +210,7 @@ public class EventNode extends HierarchicalNode
private class EventState {
private double maxDisabledTime;
private double minEnabledTime;
private boolean isInitialTrue;
private LinkedList<TriggeredEventNode> nonPersistentEvents;
public EventState()
@ -212,7 +222,14 @@ public class EventNode extends HierarchicalNode
public void reset()
{
this.maxDisabledTime = 0;
if(isInitialTrue)
{
this.maxDisabledTime = Double.NEGATIVE_INFINITY;
}
else
{
this.maxDisabledTime = 0;
}
this.minEnabledTime = Double.POSITIVE_INFINITY;
this.nonPersistentEvents.clear();
}

View file

@ -74,6 +74,7 @@ public class HierarchicalSSADirectSimulator extends HierarchicalSimulation
setCurrentTime(simProperties.getInitialTime());
ModelSetup.setupModels(this, ModelType.HSSA);
computeFixedPoint();
for(HierarchicalModel model : this.getListOfHierarchicalModels())
{
totalPropensity.addChild(model.getPropensity().getVariable());
@ -83,7 +84,7 @@ public class HierarchicalSSADirectSimulator extends HierarchicalSimulation
{
triggeredEventList =
new PriorityQueue<TriggeredEventNode>(1, new HierarchicalEventComparator());
computeEvents();
}
setupForOutput(runNumber);
@ -116,6 +117,7 @@ public class HierarchicalSSADirectSimulator extends HierarchicalSimulation
SimulationProperties simProperties = properties.getSimulationProperties();
setCurrentTime(simProperties.getInitialTime());
restoreInitialState();
computeFixedPoint();
setupForOutput(newRun);
}
@ -137,10 +139,8 @@ public class HierarchicalSSADirectSimulator extends HierarchicalSimulation
this.initialize(properties.getSimulationProperties().getRndSeed(), 1);
}
computeFixedPoint();
printTime = simProperties.getOutputStartTime();
previousTime = 0;
computeEvents();
while (currentTime.getState().getStateValue() < timeLimit && !isCancelFlag())
{
@ -153,16 +153,11 @@ public class HierarchicalSSADirectSimulator extends HierarchicalSimulation
r2 = getRandomNumberGenerator().nextDouble();
computePropensities();
totalPropensity = getTotalPropensity();
delta_t = computeNextTimeStep(r1, totalPropensity);
delta_t = Math.log(1 / r1) / totalPropensity;
nextReactionTime = currentTime + delta_t;
nextEventTime = getNextEventTime();
nextMaxTime = currentTime + maxTimeStep;
previousTime = currentTime;
if(nextMaxTime > printTime)
{
nextMaxTime = printTime;
}
if (nextReactionTime < nextEventTime && nextReactionTime < nextMaxTime)
{
@ -251,12 +246,6 @@ public class HierarchicalSSADirectSimulator extends HierarchicalSimulation
this.totalPropensity.computeFunction(0);
}
private double computeNextTimeStep(double r1, double totalPropensity)
{
return Math.log(1 / r1) / totalPropensity;
}
private void fireRateRules(double previousTime)
{

View file

@ -135,17 +135,20 @@ public class CoreSetup
{
sim.addPrintVariable(printVariable, node, modelstate.getIndex(), false);
}
node.setState(HierarchicalUtilities.createState(sim.getCollectionType(), wrapper));
modelstate.addMappingNode(compartmentID, node);
if(!compartment.isConstant())
{
modelstate.addVariable(node);
node.getState().addState(index, HierarchicalUtilities.createState(sim.getAtomicType(), wrapper));
}
else
{
node.getState().addState(index, HierarchicalUtilities.createState(StateType.SCALAR, wrapper));
}
node.setState(HierarchicalUtilities.createState(sim.getCollectionType(), wrapper));
node.getState().addState(index, HierarchicalUtilities.createState(sim.getAtomicType(), wrapper));
if (Double.isNaN(compartment.getSize()))
{
node.getState().getState(index).setInitialValue(1);
@ -176,8 +179,7 @@ public class CoreSetup
node = new VariableNode(parameterId);
modelstate.addMappingNode(parameter.getId(), node);
node.setState(HierarchicalUtilities.createState(sim.getCollectionType(), wrapper));
node.getState().addState(index, HierarchicalUtilities.createState(sim.getAtomicType(), wrapper));
String printVariable = container.getPrefix() + parameter.getId();
if(sim.getProperties().getSimulationProperties().getIntSpecies().contains(printVariable))
@ -189,17 +191,21 @@ public class CoreSetup
{
sim.addPrintVariable(printVariable, node, modelstate.getIndex(), false);
}
if(!parameter.isConstant())
{
modelstate.addVariable(node);
node.getState().addState(index, HierarchicalUtilities.createState(sim.getAtomicType(), wrapper));
}
else
{
node.getState().addState(index, HierarchicalUtilities.createState(StateType.SCALAR, wrapper));
}
if(parameter.isSetValue())
{
node.getState().getState(index).setInitialValue(parameter.getValue());
}
if(!parameter.isConstant())
{
modelstate.addVariable(node);
}
}
}
@ -231,22 +237,24 @@ public class CoreSetup
{
sim.addPrintVariable(printVariable, node, modelstate.getIndex(), isConcentration);
}
if(sim.getProperties().getSimulationProperties().getIntSpecies().size() == 0 )
else if(sim.getProperties().getSimulationProperties().getIntSpecies().size() == 0 )
{
sim.addPrintVariable(printVariable, node, modelstate.getIndex(), isConcentration);
}
if(!species.isConstant())
{
modelstate.addVariable(node);
}
boolean isBoundary = species.getBoundaryCondition();
boolean isOnlySubstance = species.getHasOnlySubstanceUnits();
node.setState(HierarchicalUtilities.createState(sim.getCollectionType(), wrapper));
node.getState().addState(index, HierarchicalUtilities.createState(sim.getAtomicType(), wrapper));
if(!species.isConstant())
{
modelstate.addVariable(node);
node.getState().addState(index, HierarchicalUtilities.createState(sim.getAtomicType(), wrapper));
}
else
{
node.getState().addState(index, HierarchicalUtilities.createState(StateType.SCALAR, wrapper));
}
node.getState().getState(index).setBoundaryCondition(isBoundary);
node.getState().getState(index).setHasOnlySubstance(isOnlySubstance);
@ -284,8 +292,16 @@ public class CoreSetup
SpeciesReferenceNode node = new SpeciesReferenceNode();
node.setState(HierarchicalUtilities.createState(sim.getCollectionType(), wrapper));
node.getState().addState(index, HierarchicalUtilities.createState(sim.getAtomicType(), wrapper));
if(!product.isConstant())
{
modelstate.addVariable(node);
node.getState().addState(index, HierarchicalUtilities.createState(sim.getAtomicType(), wrapper));
}
else
{
node.getState().addState(index, HierarchicalUtilities.createState(StateType.SCALAR, wrapper));
}
node.getState().getState(modelstate.getIndex()).setInitialValue(stoichiometryValue);
SpeciesNode species = (SpeciesNode) modelstate.getNode(product.getSpecies());
@ -323,8 +339,15 @@ public class CoreSetup
SpeciesReferenceNode node = new SpeciesReferenceNode();
node.setState(HierarchicalUtilities.createState(sim.getCollectionType(), wrapper));
node.getState().addState(index, HierarchicalUtilities.createState(sim.getAtomicType(), wrapper));
if(!reactant.isConstant())
{
modelstate.addVariable(node);
node.getState().addState(index, HierarchicalUtilities.createState(sim.getAtomicType(), wrapper));
}
else
{
node.getState().addState(index, HierarchicalUtilities.createState(StateType.SCALAR, wrapper));
}
node.getState().getState(modelstate.getIndex()).setInitialValue(stoichiometryValue);
SpeciesNode species = (SpeciesNode) modelstate.getNode(reactant.getSpecies());
@ -457,15 +480,8 @@ public class CoreSetup
boolean initValue = trigger.isSetInitialValue() ? trigger.getInitialValue() : false;
node.getState().getState(index).setUseTriggerValue(useValuesFromTrigger);
node.getState().getState(index).setPersistent(isPersistent);
node.setInitialTrue(index, initValue);
if (!initValue)
{
node.setMaxDisabledTime(index, 0);
if (node.computeTrigger(modelstate.getIndex()))
{
node.setMinEnabledTime(index, 0);
}
}
if (event.isSetPriority())
{
Priority priority = event.getPriority();