diff --git a/experiments/SIR/outbreak5.formula.xml b/experiments/SIR/outbreak5.formula.xml index 72341514..5a599b3a 100644 --- a/experiments/SIR/outbreak5.formula.xml +++ b/experiments/SIR/outbreak5.formula.xml @@ -1,14 +1,14 @@ - -1.0 -5.0 - - -I - -5.0 - + + 1.0 + 5.0 + + + I + + 5.0 + diff --git a/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/SimulatorRegistrar.java b/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/SimulatorRegistrar.java index bbe21d13..f39dc328 100644 --- a/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/SimulatorRegistrar.java +++ b/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/SimulatorRegistrar.java @@ -19,8 +19,6 @@ */ package org.sybila.parasim.computation.simulation; -import dk.ange.octave.OctaveEngine; -import dk.ange.octave.OctaveEngineFactory; import org.slf4j.Logger; import org.slf4j.LoggerFactory; import org.sybila.parasim.computation.lifecycle.api.annotations.ComputationInstanceScope; @@ -28,7 +26,6 @@ import org.sybila.parasim.computation.simulation.cpu.SimpleAdaptiveStepSimulator; import org.sybila.parasim.computation.simulation.octave.LsodeEngineFactory; import org.sybila.parasim.computation.simulation.octave.OctaveSimulationEngineFactory; -import org.sybila.parasim.computation.simulation.simulationcore.SimCoreSimulationEngineFactory; import org.sybila.parasim.core.annotation.Provide; import org.sybila.parasim.core.api.configuration.ExtensionDescriptor; import org.sybila.parasim.core.api.configuration.ExtensionDescriptorMapper; diff --git a/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/cpu/SimulationEngine.java b/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/cpu/SimulationEngine.java index 0f33bbc7..1a98035a 100644 --- a/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/cpu/SimulationEngine.java +++ b/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/cpu/SimulationEngine.java @@ -29,14 +29,18 @@ */ public interface SimulationEngine { + /** + * Closes SimulationEngine. This method is called after the simulation of the whole space is performed. + * Enables to close external software e.g. GNU Octave + */ void close(); /** - * Performs simulation of a differential equation system from given point in time to timelimit - * @param point where to start simulation, contains start time + * Performs simulation of a differential equation system from given point. + * @param point initial point where to start simulation, contains start time and values of state variables and parameters * @param odeSystem differential equation system * @param timeLimit end time of simulation - * @param configuration relative error, absolute error and time step configuration + * @param configuration contains relative error, absolute error and time step configuration * @return computed trajectory */ Trajectory simulate(Point point, OdeSystem odeSystem, double timeLimit, PrecisionConfiguration configuration); diff --git a/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/cpu/SimulationEngineFactory.java b/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/cpu/SimulationEngineFactory.java index 07e45a13..0098b6e7 100644 --- a/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/cpu/SimulationEngineFactory.java +++ b/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/cpu/SimulationEngineFactory.java @@ -28,8 +28,17 @@ public interface SimulationEngineFactory{ boolean isAvailable(); + /** + * To get simulation engine for the thread which this method is called from. + * New SimulationeEngine instance is constructed, if none exists according to THREAD_SIMULATION_ENGINE_MAP. + * @param stepLimit max number of simulation steps of the simulation engine + * @return SimulationEngine + */ E simulationEngine(long stepLimit); + /** + * Hash map that allows closing of all SimulationEngine instances after the simulation of the whole space finishes + */ ConcurrentHashMap THREAD_SIMULATION_ENGINE_MAP = new ConcurrentHashMap<>(); } diff --git a/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/octave/LsodeEngineFactory.java b/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/octave/LsodeEngineFactory.java index f4797828..d0881438 100644 --- a/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/octave/LsodeEngineFactory.java +++ b/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/octave/LsodeEngineFactory.java @@ -23,10 +23,8 @@ import dk.ange.octave.OctaveEngineFactory; import dk.ange.octave.type.OctaveDouble; import java.util.Arrays; -import java.util.function.Function; import org.sybila.parasim.computation.simulation.api.PrecisionConfiguration; -import org.sybila.parasim.computation.simulation.cpu.SimulationEngine; import org.sybila.parasim.computation.simulation.cpu.SimulationEngineFactory; import org.sybila.parasim.model.ode.OctaveOdeSystem; import org.sybila.parasim.model.trajectory.Point; @@ -47,7 +45,7 @@ public LsodeEngineFactory(IntegrationMethod integrationMethod) { public boolean isAvailable() { try { new OctaveEngineFactory().getScriptEngine(); -// octave.close(); //closing octave after opening? but this method seems to work not +// octave.close(); //should close octave after opening? -- luckily this method is not called anywhere (old one) return true; } catch (Exception ignored) { return false; @@ -55,7 +53,7 @@ public boolean isAvailable() { } @Override - public OctaveSimulationEngine simulationEngine(long stepLimit) {//TODO correct retyping? + public OctaveSimulationEngine simulationEngine(long stepLimit) { return (OctaveSimulationEngine) SimulationEngineFactory.THREAD_SIMULATION_ENGINE_MAP.computeIfAbsent(Thread.currentThread(), thread -> new LsodeEngine(new OctaveEngineFactory().getScriptEngine(), integrationMethod, stepLimit)); } diff --git a/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/octave/OctaveSimulationEngine.java b/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/octave/OctaveSimulationEngine.java index 9f981202..bd7217dc 100644 --- a/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/octave/OctaveSimulationEngine.java +++ b/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/octave/OctaveSimulationEngine.java @@ -23,13 +23,10 @@ import dk.ange.octave.type.OctaveDouble; import java.util.ArrayList; import java.util.List; -import java.util.concurrent.atomic.AtomicLong; - import org.apache.commons.lang3.Validate; import org.slf4j.Logger; import org.slf4j.LoggerFactory; import org.sybila.parasim.computation.simulation.api.PrecisionConfiguration; -import org.sybila.parasim.computation.simulation.cpu.SimpleAdaptiveStepSimulator; import org.sybila.parasim.computation.simulation.cpu.SimulationEngine; import org.sybila.parasim.model.math.Parameter; import org.sybila.parasim.model.math.ParameterValue; @@ -40,8 +37,6 @@ import org.sybila.parasim.model.trajectory.Point; import org.sybila.parasim.model.trajectory.Trajectory; -import javax.management.ServiceNotFoundException; - /** * @author Jan Papousek */ @@ -63,45 +58,27 @@ public OctaveSimulationEngine(OctaveEngine octave, long stepLimit) { /** * caching simulation time (can save time if the last simulation lasts the same time as the previous) */ - float[] l_times; + private float[] l_times; /** * start time of last simulation */ - float l_startTime; + private float l_startTime; /** * end of last simulation */ - double l_endTime; + private double l_endTime; /** * time step of last simulation */ - float timeStep; + private float timeStep; @Override public void close() { getOctave().close(); } - private static AtomicLong totalSettingTime; - private static AtomicLong totalSimulationTime; - private static AtomicLong totalParsingTime; - - private static AtomicLong numOfSimulations; - - static { - totalSettingTime = new AtomicLong(); - totalSettingTime.set(0); - totalSimulationTime = new AtomicLong(); - totalSimulationTime.set(0); - totalParsingTime = new AtomicLong(); - totalParsingTime.set(0); - numOfSimulations = new AtomicLong(); - numOfSimulations.set(0); - } - @Override public Trajectory simulate(Point point, OdeSystem odeSystem, double timeLimit, PrecisionConfiguration precision) { - long settingStartTime = System.nanoTime(); // load parameter values List paramValues = loadParameterValues(point, odeSystem); // create substituted octave ode system @@ -111,24 +88,13 @@ public Trajectory simulate(Point point, OdeSystem odeSystem, double timeLimit, P if (numOfIterations > getStepLimit()) { throw new IllegalStateException("Can't simulate the trajectory because the number of iterations <" + numOfIterations + "> is higher than the given limit <" + getStepLimit() + ">."); } - long settingTime = System.nanoTime() - settingStartTime; - - long simulationStartTime = System.nanoTime(); double[] loadedData = rawSimulation(point, octaveOdeSystem, numOfIterations, precision).getData(); - long simulationTime = System.nanoTime() - simulationStartTime; - - - long parsingStartTime = System.nanoTime(); float[] data = new float[loadedData.length]; for (int dim = 0; dim < octaveOdeSystem.dimension(); dim++) { for (int i = 0; i < loadedData.length / octaveOdeSystem.dimension(); i++) { data[dim + i * octaveOdeSystem.dimension()] = (float) loadedData[dim * (loadedData.length / octaveOdeSystem.dimension()) + i]; } } - long parsingTime = System.nanoTime() - parsingStartTime; - - long timeStartTime = System.nanoTime(); - //if simulation time changed, rewrite cached time array if(l_startTime != point.getTime() || l_endTime != timeLimit || timeStep != precision.getTimeStep()){ l_startTime = point.getTime(); @@ -141,23 +107,6 @@ public Trajectory simulate(Point point, OdeSystem odeSystem, double timeLimit, P l_times[i] = time; } } - - long timeTime = System.nanoTime() - timeStartTime; - - totalSettingTime.addAndGet(settingTime); - totalSimulationTime.addAndGet(simulationTime); - totalParsingTime.addAndGet(parsingTime + timeTime); - numOfSimulations.addAndGet(1); - - System.out.println("AVERAGE TIME"); - System.out.printf("Setting time: %.9f ms\n", (totalSettingTime.get()/numOfSimulations.get()) / 1000000000.0); - System.out.printf("Simulation time: %.9f ms\n", (totalSimulationTime.get()/numOfSimulations.get()) / 1000000000.0); -// System.out.printf("Time time: %.9f ms\n", timeTime / 1000000000.0); - System.out.printf("Parsing time: %.9f ms\n", (totalParsingTime.get()/numOfSimulations.get()) / 1000000000.0); - System.out.println("Number of simulated trajectories: " + numOfSimulations.get()); - - - if (paramValues.isEmpty()) { return new ArrayTrajectory(data, l_times, point.getDimension()); } else { diff --git a/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/simulationcore/SimCoreSimulationEngine.java b/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/simulationcore/SimCoreSimulationEngine.java index 074c2866..40bfc599 100644 --- a/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/simulationcore/SimCoreSimulationEngine.java +++ b/extensions/computation-simulation-impl/src/main/java/org/sybila/parasim/computation/simulation/simulationcore/SimCoreSimulationEngine.java @@ -3,22 +3,17 @@ import org.apache.commons.math.ode.DerivativeException; import org.sbml.jsbml.Model; import org.sbml.jsbml.validator.ModelOverdeterminedException; -import org.simulator.math.odes.AbstractDESSolver; import org.simulator.math.odes.AdaptiveStepsizeIntegrator; import org.simulator.math.odes.MultiTable; import org.simulator.math.odes.RosenbrockSolver; import org.simulator.sbml.SBMLinterpreter; import org.sybila.parasim.computation.simulation.api.PrecisionConfiguration; -import org.sybila.parasim.computation.simulation.api.SimulationException; import org.sybila.parasim.computation.simulation.cpu.SimulationEngine; import org.sybila.parasim.model.math.Parameter; import org.sybila.parasim.model.math.Variable; import org.sybila.parasim.model.ode.OdeSystem; import org.sybila.parasim.model.trajectory.*; -import java.util.concurrent.atomic.AtomicLong; - - /** * @author Vojtech Bruza */ @@ -31,51 +26,27 @@ public void close() { /** * array for caching time of simulations for Parasim */ - float[] l_times; + private float[] l_times; /** * array for caching time of simulations for Simulatin Core Library */ - double [] l_timesDouble; + private double [] l_timesDouble; /** * start of last simulation */ - float l_startTime; + private float l_startTime; /** * end of last simulation */ - double l_endTime; + private double l_endTime; /** * time step of previous simulation */ - float timeStep; - - - private static AtomicLong totalSettingTime; - private static AtomicLong totalSimulationTime; - private static AtomicLong totalParsingTime; - - private static AtomicLong numOfSimulations; - - static { - totalSettingTime = new AtomicLong(); - totalSettingTime.set(0); - totalSimulationTime = new AtomicLong(); - totalSimulationTime.set(0); - totalParsingTime = new AtomicLong(); - totalParsingTime.set(0); - numOfSimulations = new AtomicLong(); - numOfSimulations.set(0); - } - + private float timeStep; @Override public Trajectory simulate(Point point, OdeSystem odeSystem, double timeLimit, PrecisionConfiguration precision) { - long settingStartTime = System.nanoTime(); Model model = odeSystem.getOriginalModel(); - -// System.out.println("PARAMETER VALUES: " + odeSystem.getAvailableParameters().keySet() +" VARIABLES: " + odeSystem.getVariables().keySet()); -// // DONE how to recognize if the index is same as in my Model (corresponding parameters and variables) - // //SETTING VARIABLES for(Variable variable : odeSystem.getVariables().values()){ if (!variable.isSubstituted()) { @@ -83,37 +54,19 @@ public Trajectory simulate(Point point, OdeSystem odeSystem, double timeLimit, P model.getSpecies(variable.getName()).setValue(point.getValue(variable.getIndex())); } } - -// //odeSystem.getAvailableParameters() returns only parameters, use odeSystem.getVariables() to get variables -// //DONE findOut if parameters in parasim are also variables (if it is synonym) - -// //SETTING PARAMETERS -// for(Parameter parameter : odeSystem.getAvailableParameters().values()){ -// if (!parameter.isSubstituted()) { -// System.out.println(parameter.getName() + " - VALUE: " + point.getValue(parameter.getIndex())); -// //set parameters values in model -// model.getParameter(parameter.getName()).setValue(point.getValue(parameter.getIndex()));//what is the difference?? "parameter.evaluate(point)" - doesnt work if parameter is not substituted vs "point.getValue(parameter.getIndex())" - works -// } -// } -// //DONE set parameter value in model according to ode system parameter value -// //DONE create substituted odeSystem (model) - - //DONE!! SET PARAMETERS + //SETTING PARAMETERS for(Parameter parameter : odeSystem.getAvailableParameters().values()){ if (!parameter.isSubstituted()) { //substituted are those that are set at the beginning (not perturbating over them) //set parameters values in model model.getParameter(parameter.getName()).setValue(point.getValue(parameter.getIndex()));//what is the difference?? "parameter.evaluate(point)" - doesn't work if parameter is not substituted vs "point.getValue(parameter.getIndex())" - works } } - - //SET SIMULATION TIME (NUM OF ITERATIONS) - START, END, NUM OF ITERATION, TIME STEP long numOfIterations = Math.round(Math.ceil((1.05 * timeLimit - point.getTime()) / precision.getTimeStep())); //magical constant 1.05 taken from Jan Papousek's implementation => guessing it is for rendering reasons (to simulate a bit more data than a user expects) if (numOfIterations > Integer.MAX_VALUE) { // step limit is max integer value, because it needs to be retyped from long later on throw new IllegalStateException("Can't simulate the trajectory because the number of iterations <" + numOfIterations + "> has exceeded its maximum <" + Integer.MAX_VALUE + ">."); } - - //TIME - needs to be float due to communication with Parasim and double due to SimCore + //TIME - needs to be float due to communication with other modules of Parasim and double due to SimCore if(l_startTime != point.getTime() || l_endTime != timeLimit || timeStep != precision.getTimeStep()){ l_startTime = point.getTime(); l_endTime = timeLimit; @@ -127,8 +80,7 @@ public Trajectory simulate(Point point, OdeSystem odeSystem, double timeLimit, P l_timesDouble[j] = timeFloat; } } - - //SIMULATION METHOD CHOICE + //SET SIMULATION METHOD SBMLinterpreter interpreter = null; try { interpreter = new SBMLinterpreter(model); @@ -148,35 +100,12 @@ public Trajectory simulate(Point point, OdeSystem odeSystem, double timeLimit, P } } solver.setAbsTol(maxAbsoluteError); - long settingTime = System.nanoTime() - settingStartTime; - //SIMULATION - long simulationStartTime = System.nanoTime(); try { solution = solver.solve(interpreter, interpreter.getInitialValues(), l_timesDouble); - } catch (DerivativeException e) { + } catch (DerivativeException | NullPointerException e) { e.printStackTrace(); } - long simulationTime = System.nanoTime() - simulationStartTime; - //PARSING DATA TO TRAJECTORY - - //help printing -// for (int i = 0; i < solution.getBlockCount(); i++) { -// System.out.println("block" + i + " name: " + solution.getBlock(i).getName()); -// for (int c = 1; c < solution.getBlock(i).getColumnCount(); c++) { -// System.out.println("\tcolumn" + c + " name: " + solution.getBlock(i).getColumn(c).getColumnName()); -// } -// } -// System.out.println("COLUMNS:"); -// for (Parameter p : odeSystem.getAvailableParameters().values()){ -// System.out.println("parameter: " + p.getName() + " column: " + solution.getColumn(p.getName()).getColumnName()); -// } -// for (Variable v : odeSystem.getVariables().values()){ -// System.out.println("variable: " + v.getName() + " column: " + solution.getColumn(v.getName()).getColumnName()); -// } -// System.out.println("Start time: " + point.getTime() + " End Time: " + timeLimit); -// System.out.println("Number of iterations: " + numOfIterations); - long parsingStartTime = System.nanoTime(); int numberOfSteps = 0; try{ numberOfSteps = solution.getRowCount(); @@ -189,20 +118,6 @@ public Trajectory simulate(Point point, OdeSystem odeSystem, double timeLimit, P simulatedData[currentVariable + i * odeSystem.getVariables().size()] = (float) solution.getColumn(odeSystem.getVariable(currentVariable).getName()).getValue(i); } } - long parsingTime = System.nanoTime() - parsingStartTime; - - totalSettingTime.addAndGet(settingTime); - totalSimulationTime.addAndGet(simulationTime); - totalParsingTime.addAndGet(parsingTime); - numOfSimulations.addAndGet(1); - - System.out.println("AVERAGE TIME"); - System.out.printf("Setting time: %.9f ms\n", (totalSettingTime.get()/numOfSimulations.get()) / 1000000000.0); - System.out.printf("Simulation time: %.9f ms\n", (totalSimulationTime.get()/numOfSimulations.get()) / 1000000000.0); -// System.out.printf("Time time: %.9f ms\n", timeTime / 1000000000.0); - System.out.printf("Parsing time: %.9f ms\n", (totalParsingTime.get()/numOfSimulations.get()) / 1000000000.0); - System.out.println("Number of simulated trajectories: " + numOfSimulations.get()); - //OUTPUT TRAJECTORY if (odeSystem.getAvailableParameters().isEmpty()) { return new ArrayTrajectory(simulatedData, l_times, point.getDimension());