Merge pull request #178 from StochSS/tspan-fix

Timespan Fix

Author: briandrawert

spatialpy/Model.py

@@ -95,7 +95,7 @@ def __init__(self, name="spatialpy"):
 
         ######################
         self.tspan = None
-        self.timestep_size = 1e-5
+        self.timestep_size = None
         self.num_timesteps = None
         self.output_freq = None
 
@@ -188,51 +188,50 @@ def __check_if_complete(self):
 
     def set_timesteps(self, output_interval, num_steps, timestep_size=None):
         """" Set the simlation time span parameters
-        Args:
-            :param output_interval: size of each output timestep in seconds
-            :type output_interval:  float
-            :param num_steps: total number of steps to take. Note: the number of output times will be num_steps+1 as the first
-              output will be at time zero.
-            :type num_steps: int
-            :param step_size: Size of each timestep in seconds
-            :type step_size: float
-            :param num_steps: Total number of steps to take
-            :type num_steps: int
-
+        :param output_interval: size of each output timestep in seconds
+        :type output_interval:  float
+        :param num_steps: total number of steps to take. Note: the number of output times will be num_steps+1 as the first
+          output will be at time zero.
+        :type num_steps: int
+        :param timestep_size: Size of each timestep in seconds
+        :type timestep_size: float
         """
         if timestep_size is not None:
             self.timestep_size = timestep_size
         if self.timestep_size is None:
-            raise ModelError("timestep_size is not set")
-
-        self.output_freq = math.ceil(output_interval/self.timestep_size)
+            self.timestep_size = output_interval
+        
+        self.output_freq = output_interval/self.timestep_size
+        if self.output_freq < self.timestep_size:
+            raise ModelError("Timestep size exceeds output frequency.")
 
         self.num_timesteps = math.ceil(num_steps * self.output_freq)
-        # array of step numbers corresponding to the simulation times in the timespan
-        self.timespan_steps = numpy.linspace(0,self.num_timesteps,
-                        num=math.ceil(self.num_timesteps/self.output_freq)+1, dtype=int)
-
-        self.tspan = numpy.linspace(0,self.num_timesteps*self.timestep_size,self.num_timesteps)
 
+        # array of step numbers corresponding to the simulation times in the timespan
+        output_steps = numpy.arange(0, self.num_timesteps + self.timestep_size, self.output_freq)
+        self.output_steps = numpy.unique(numpy.round(output_steps).astype(int))
+        sim_steps = numpy.arange(0, self.num_timesteps + self.timestep_size, self.timestep_size)
+        self.tspan = numpy.zeros((self.output_steps.size), dtype=float)
+        for i, step in enumerate(self.output_steps):
+            self.tspan[i] = sim_steps[step]
+    
     def timespan(self, time_span, timestep_size=None):
         """
         Set the time span of simulation. The SSA-SDPD engine does not support
         non-uniform timespans.
 
         :param tspan: Evenly-spaced list of times at which to sample the species populations during the simulation.
         :type tspan: numpy.ndarray
+        :param timestep_size: Size of each timestep in seconds
+        :type timestep_size: float
         """
 
-        self.tspan = time_span
-        if timestep_size is not None:
-            self.timestep_size = timestep_size
-
         items_diff = numpy.diff(time_span)
         items = map(lambda x: round(x, 10), items_diff)
         isuniform = (len(set(items)) == 1)
 
         if isuniform:
-            self.set_timesteps( items_diff[0], len(items_diff) )
+            self.set_timesteps(items_diff[0], len(items_diff), timestep_size=timestep_size)
         else:
             raise ModelError("Only uniform timespans are supported")
 

spatialpy/Result.py

@@ -193,10 +193,7 @@ def read_step(self, step_num, debug=False):
 
         if debug:
             print("read_step({0}) ".format(step_num), end='')
-        #num = int(step_num * self.model.output_freq)
-        #num = int(step_num / self.model.timestep_size)
-        num = self.model.timespan_steps[step_num]
-        filename = os.path.join(self.result_dir, "output{0}.vtk".format(num))
+        filename = os.path.join(self.result_dir, "output{0}.vtk".format(step_num))
 
         if debug:
             print("opening '{0}'".format(filename))
@@ -238,8 +235,7 @@ def get_timespan(self):
 
         """
 
-        self.tspan = numpy.linspace(0,self.model.num_timesteps,
-                num=math.ceil(self.model.num_timesteps/self.model.output_freq)+1) * self.model.timestep_size
+        self.tspan = self.model.tspan
         return self.tspan
 
     def get_species(self, species, timepoints=None, concentration=False, deterministic=False, debug=False):
@@ -278,13 +274,9 @@ def get_species(self, species, timepoints=None, concentration=False, determinist
         if spec_name not in self.model.listOfSpecies.keys():
             raise ResultError("Species '{0}' not found".format(spec_name))
 
-        #t_index_arr = numpy.linspace(0,self.model.num_timesteps,
-        #                    num=self.model.num_timesteps+1, dtype=int)
-        #t_index_arr = numpy.linspace(0,self.model.num_timesteps,
-        #        num=math.ceil(self.model.num_timesteps/self.model.output_freq)+1) 
         lt=len(self.get_timespan())-1
         t_index_arr = numpy.linspace(0,lt,num=lt+1,dtype=int)
-
+        
         if timepoints is not None:
             if isinstance(timepoints,float):
                 raise ResultError("timepoints argument must be an integer, the index of time timespan")
@@ -572,8 +564,9 @@ def get_property(self, property_name, timepoints=None):
             :rtype: numpy.ndarray
         """
 
-        t_index_arr = numpy.linspace(0,self.model.num_timesteps,
-                            num=self.model.num_timesteps+1, dtype=int)
+        lt=len(self.get_timespan())-1
+        t_index_arr = numpy.linspace(0,lt,num=lt+1,dtype=int)
+        
         num_voxel = self.model.domain.get_num_voxels()
 
         if timepoints is not None:

spatialpy/Solver.py

@@ -119,6 +119,11 @@ def compile(self, debug=False, profile=False):
             :type profile: bool
         """
 
+        if self.model.timestep_size is None:
+            self.model.timestep_size = 1e-5
+        if self.model.output_freq < self.model.timestep_size:
+            raise ModelError("Timestep size exceeds output frequency.")
+
         # Create the models expression utility
         self.model.get_expression_utility()
 
@@ -580,7 +585,6 @@ def __create_propensity_file(self, file_name=None):
 
         system_config += "system->dt = {0};\n".format(self.model.timestep_size)
         system_config += "system->nt = {0};\n".format(self.model.num_timesteps)
-        system_config += "system->output_freq = {0};\n".format(self.model.output_freq)
         if self.h is None:
             self.h = self.model.domain.find_h()
         if self.h == 0.0:
@@ -625,6 +629,15 @@ def __create_propensity_file(self, file_name=None):
             init_bc += bc.expression()
         propfilestr = propfilestr.replace("__BOUNDARY_CONDITIONS__", init_bc)
 
+        output_step = "unsigned int get_next_output(ParticleSystem* system)\n{\n"
+        output_step += "static int index = 0;\n"
+        output_step += "const std::vector<unsigned int> output_steps = {"
+        output_step += f"{', '.join(self.model.output_steps.astype(str).tolist())}"
+        output_step += "};\nunsigned int next_step = output_steps[index];\n"
+        output_step += "index++;\n"
+        output_step += "return next_step;\n}\n"
+        propfilestr = propfilestr.replace("__DEFINE_GET_NEXT_OUTPUT__", output_step)
+
         #### Write the data to the file ####
         propfile.write(propfilestr)
         propfile.close()

spatialpy/ssa_sdpd-c-simulation-engine/include/particle_system.hpp

@@ -61,7 +61,6 @@ namespace Spatialpy{
         double dt;
         unsigned int nt;
         unsigned int current_step;
-        unsigned int output_freq;
         double xlo, xhi, ylo, yhi, zlo, zhi;
         double h;
         double c0;

spatialpy/ssa_sdpd-c-simulation-engine/include/propensities.hpp

@@ -36,6 +36,7 @@ namespace Spatialpy{
     /* Global variable that can be used to pass parameters to the propensity functions. */
     extern double *parameters;
 
+    unsigned int get_next_output(ParticleSystem* system);
     /* Definition of the propensity function. */
     // double rfun(const int *x, double t, const double vol, const double *data, int sd, int voxel, int *xx, const size_t *irK, const size_t *jcK, const double *prK)
     //typedef double (*PropensityFun)(const int *, double, double, const double *, int, int, int *, const size_t *, const size_t *, const double *);

spatialpy/ssa_sdpd-c-simulation-engine/propensity_file_template.cpp

@@ -41,6 +41,8 @@ namespace Spatialpy{
     __DEFINE_REACTIONS__
     /* Deterministic RHS definitions */
     __DEFINE_CHEM_FUNS__
+    /* Definition for get next output step */
+    __DEFINE_GET_NEXT_OUTPUT__
 
     PropensityFun *ALLOC_propensities(void)
     {

spatialpy/ssa_sdpd-c-simulation-engine/src/output.c

@@ -106,6 +106,7 @@ namespace Spatialpy{
         int i;
         char filename[256];
         int np = output_buffer_current_num_particles;
+        static unsigned int output_index = 0;
         if(output_buffer_current_step == 0){
             sprintf(filename, "output0_boundingBox.vtk");
             if(debug_flag){printf("Writing file '%s'\n", filename);}
@@ -125,7 +126,7 @@ namespace Spatialpy{
             fprintf(fp, "%lf %lf\n", system->zlo, system->zhi);
             fclose(fp);
         }
-        sprintf(filename,"output%u.vtk",output_buffer_current_step);
+        sprintf(filename,"output%u.vtk", output_index++);
         if(debug_flag){printf("Writing file '%s'\n", filename);}
         if((fp = fopen(filename,"w+"))==NULL){
             perror("Can't write output vtk file");exit(1);

spatialpy/ssa_sdpd-c-simulation-engine/src/simulate_threads.c

@@ -239,7 +239,7 @@ namespace Spatialpy{
                 // Release output thread
                 if(debug_flag) printf("[%i] Starting the Output threads\n",system->current_step);
                 pthread_barrier_wait(&begin_output_barrier);
-                next_output_step += system->output_freq;
+                next_output_step = get_next_output(system);
                 // Wait until output threads are done
                 pthread_barrier_wait(&end_output_barrier);
                 if(debug_flag) printf("[%i] Output threads finished\n",system->current_step);