Merge pull request #288 from StochSS/staging

Staging

Author: BryanRumsey

.github/workflows/pypi_publish.yaml

@@ -0,0 +1,30 @@
+name: Publish SpatialPy
+
+on:
+  release:
+    types: [published]
+
+jobs:
+  deploy:
+    runs-on: ubuntu-latest
+
+    steps:
+    - uses: actions/checkout@v2
+    - name: Set up Python
+      uses: actions/setup-python@v2
+      with:
+        python-version: '3.7'
+        
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install build
+        
+    - name: Build package
+      run: python -m build
+      
+    - name: Publish package
+      uses: pypa/gh-action-pypi-publish@27b31702a0e7fc50959f5ad993c78deac1bdfc29
+      with:
+        user: __token__
+        password: ${{ secrets.PYPI_API_TOKEN }}

.graphics/birth-death-example-plot.png

@@ -74,6 +74,68 @@ SpatialPy provides simple object-oriented abstractions for defining a model of a
 
 The `run()` method can be customized using keyword arguments to select different solvers, random seed, data return type and more.  For more detailed examples on how to use SpatialPy, please see the Jupyter notebooks contained in the [examples](https://github.com/StochSS/SpatialPy/tree/main/examples) subdirectory.
 
+### _Simple example to illustrate the use of SpatialPy_
+
+In SpatialPy, a model is expressed as an object. Components, such as the domains, reactions, biochemical species, and characteristics such as the time span for simulation, are all defined within the model. The following Python code represents our spatial birth death model using SpatialPy's facility:
+
+```python
+def create_birth_death(parameter_values=None):
+    # First call the gillespy2.Model initializer.
+    model = spatialpy.Model(name='Spatial Birth-Death')
+    
+    # Define Domain Type IDs as constants of the Model
+    model.HABITAT = "Habitat"
+
+    # Define domain points and attributes of a regional space for simulation.
+    domain = spatialpy.Domain.create_2D_domain(
+        xlim=(0, 1), ylim=(0, 1), nx=10, ny=10, type_id=model.HABITAT, fixed=True
+    )
+    model.add_domain(domain)
+
+    # Define variables for the biochemical species representing Rabbits.
+    Rabbits = spatialpy.Species(name='Rabbits', diffusion_coefficient=0.1)
+    model.add_species(Rabbits)
+
+    # Scatter the initial condition for Rabbits randomly over all types.
+    init_rabbit_pop = spatialpy.ScatterInitialCondition(species='Rabbits', count=100)
+    model.add_initial_condition(init_rabbit_pop)
+
+    # Define parameters for the rates of creation and destruction.
+    kb = spatialpy.Parameter(name='k_birth', expression=10)
+    kd = spatialpy.Parameter(name='k_death', expression=0.1)
+    model.add_parameter([kb, kd])
+
+    # Define reactions channels which cause the system to change over time.
+    # The list of reactants and products for a Reaction object are each a
+    # Python dictionary in which the dictionary keys are Species objects
+    # and the values are stoichiometries of the species in the reaction.
+    birth = spatialpy.Reaction(name='birth', reactants={}, products={"Rabbits":1}, rate="k_birth")
+    death = spatialpy.Reaction(name='death', reactants={"Rabbits":1}, products={}, rate="k_death")
+    model.add_reaction([birth, death])
+
+    # Set the timespan of the simulation.
+    tspan = spatialpy.TimeSpan.linspace(t=10, num_points=11, timestep_size=1)
+    model.timespan(tspan)
+    return model
+```
+
+Given the model creation function above, the model can be simulated by first instantiating the model object, and then invoking the run() method on the object. The following code will run the model once to produce a sample trajectory:
+
+```python
+model = create_birth_death()
+results = model.run()
+```
+
+The results are then stored in a class `Results` object for single trajectory or for multiple trajectories.  Results can be plotted with plotly (offline) using `plot_species()` or in matplotlib using `plot_species(use_matplotlib=True)`.  For additional plotting options such as plotting from a selection of species, or statistical plotting, please see the documentation.:
+
+```python
+results.plot_species(species='Rabbits', t_val=10, use_matplotlib=True)
+```
+
+<p align="center">
+<img width="500px" src="https://raw.githubusercontent.com/StochSS/SpatialPy/readme-example/.graphics/birth-death-example-plot.png">
+</p>
+
 <!-- 
 ### Docker environment (DOES NOT WORK)
 

README.md

@@ -0,0 +1,12 @@
+Examples of using SpatialPy
+===========================
+
+The files in this directory are runnable Python examples of using [SpatialPy](https://github.com/StochSS/SpatialPy) to perform spatial deterministic/stochastic reaction-diffusion-advection simulations.  In terms of biology, they are overly simplistic and do not capture the real-life complexity of the process being modeled – the aim is not biological realism but rather to illustrate basic usage of [SpatialPy](https://github.com/StochSS/SpatialPy).
+
+* [Start Here](Start_Here.ipynb) – a [Jupyter Notebook](https://jupyter-notebook.readthedocs.io/en/stable/) demonstrating the use of SpatialPy on a simple spatial Birth Death model.
+
+* [3D Cylinder Demo](3D_Cylinder_Demo.ipynb) – a [Jupyter Notebook](https://jupyter-notebook.readthedocs.io/en/stable/) demonstrating reaction diffusion simulations using SpatialPy.
+
+* [Gravity](Gravity.ipynb) and [Weir](Weir.ipynb) – are [Jupyter Notebook](https://jupyter-notebook.readthedocs.io/en/stable/) demonstrating fluid dynamics simulations using SpatialPy.
+
+Full documentation for SpatialPy can be found at https://stochss.github.io/SpatialPy/index.html

examples/README.md

@@ -21,7 +21,7 @@
 # @website https://github.com/StochSS/SpatialPy
 # =============================================================================
 
-__version__      = '1.0.3'
+__version__      = '1.0.4'
 __title__        = 'SpatialPy'
 __description__  = 'Python Interface for Spatial Stochastic Biochemical Simulations'
 __url__          = 'https://spatialpy.github.io/SpatialPy/'

spatialpy/__version__.py

@@ -87,6 +87,8 @@ def __init__(self, xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=N
         if type_id is not None and not isinstance(type_id, (int, str)):
             raise BoundaryConditionError("Type-ID must be of type int.")
         elif type_id is not None:
+            if "UnAssigned" in type_id:
+                raise BoundaryConditionError("'UnAssigned' is not a valid type_id")
             type_id = f"type_{type_id}"
         if target is None or not (isinstance(target, (str, Species)) or
                 type(target).__name__ == 'Species' or property in ('nu', 'rho', 'v')):

spatialpy/core/boundarycondition.py

@@ -68,13 +68,13 @@ def __init__(self, numpoints, xlim, ylim, zlim, rho0=1.0, c0=10, P0=None, gravit
 
         self.vol = numpy.zeros((numpoints), dtype=float)
         self.mass = numpy.zeros((numpoints), dtype=float)
-        self.type_id = numpy.array([None] * numpoints, dtype=object)
+        self.type_id = numpy.array(["type_UnAssigned"] * numpoints, dtype=object)
         self.nu = numpy.zeros((numpoints), dtype=float)
         self.c = numpy.zeros((numpoints), dtype=float)
         self.rho = numpy.zeros((numpoints), dtype=float)
         self.fixed = numpy.zeros((numpoints), dtype=bool)
         self.listOfTypeIDs = []
-        self.typeNdxMapping = OrderedDict()
+        self.typeNdxMapping = OrderedDict({"type_UnAssigned": 0})
         self.typeNameMapping = None
 
         self.rho0 = rho0
@@ -136,7 +136,7 @@ def compile_prep(self):
 
         :raises DomainError: If a type_id is not set or rho=0 for a particle.
         """
-        if self.type_id.tolist().count(None) > 0:
+        if self.type_id.tolist().count("type_UnAssigned") > 0:
             raise DomainError(f"Particles must be assigned a type_id.")
         if numpy.count_nonzero(self.rho) < len(self.rho):
             raise DomainError(f"Rho must be a positive value.")
@@ -184,7 +184,10 @@ def add_point(self, point, vol=0, mass=0, type_id=1, nu=0, fixed=False, rho=None
             if (char in string.punctuation and char != "_") or char == " ":
                 raise DomainError(f"Type_id cannot contain {char}")
         if type_id not in self.typeNdxMapping:
-            self.typeNdxMapping[type_id] = len(self.typeNdxMapping) + 1
+            if "UnAssigned" in type_id:
+                self.typeNdxMapping[type_id] = 0
+            else:
+                self.typeNdxMapping[type_id] = len(self.typeNdxMapping)
 
         if rho is None:
             rho = mass / vol
@@ -240,7 +243,10 @@ def set_properties(self, geometry_ivar, type_id, vol=None, mass=None, nu=None, r
             if (char in string.punctuation and char != "_") or char == " ":
                 raise DomainError(f"Type_id cannot contain '{char}'")
         if type_id not in self.typeNdxMapping:
-            self.typeNdxMapping[type_id] = len(self.typeNdxMapping) + 1
+            if "UnAssigned" in type_id:
+                self.typeNdxMapping[type_id] = 0
+            else:
+                self.typeNdxMapping[type_id] = len(self.typeNdxMapping)
         # apply the type to all points, set type for any points that match
         count = 0
         on_boundary = self.find_boundary_points()
@@ -571,6 +577,9 @@ def plot_types(self, width=None, height=None, colormap=None, size=None, title=No
         :returns: Plotly figure of domain types if, use_matplotlib=False and return_plotly_figure=True
         :rtype: None or dict
         '''
+        if len(self.vertices) == 0:
+            raise DomainError("The domain does not contain particles.")
+
         from spatialpy.core.result import _plotly_iterate # pylint: disable=import-outside-toplevel
 
         if not use_matplotlib:
@@ -807,7 +816,10 @@ def read_stochss_subdomain_file(self, filename, type_ids=None):
                         if (char in string.punctuation and char != "_") or char == " ":
                             raise DomainError(f"Type_id cannot contain {char}")
                     if type_id not in self.typeNdxMapping:
-                        self.typeNdxMapping[type_id] = len(self.typeNdxMapping) + 1
+                        if "UnAssigned" in type_id:
+                            self.typeNdxMapping[type_id] = 0
+                        else:
+                            self.typeNdxMapping[type_id] = len(self.typeNdxMapping)
 
                     self.type_id[int(ndx)] = type_id
 
@@ -834,13 +846,15 @@ def read_stochss_domain(cls, filename):
             obj = Domain(0, tuple(domain['x_lim']), tuple(domain['y_lim']), tuple(domain['z_lim']),
                         rho0=domain['rho_0'], c0=domain['c_0'], P0=domain['p_0'], gravity=domain['gravity'])
 
-            for particle in domain['particles']:
+            for i, particle in enumerate(domain['particles']):
                 try:
                     type_id = list(filter(
                         lambda d_type, t_ndx=particle['type']: d_type['typeID'] == t_ndx, domain['types']
                     ))[0]['name']
                 except IndexError:
                     type_id = particle['type']
+                if type_id == "Un-Assigned" or type_id == 0:
+                    type_id = "UnAssigned"
                 # StochSS backward compatability check for rho
                 rho = None if "rho" not in particle.keys() else particle['rho']
                 # StochSS backward compatability check for c
@@ -906,7 +920,7 @@ def create_3D_domain(cls, xlim, ylim, zlim, nx, ny, nz, type_id=1, mass=1.0,
         x_list = numpy.linspace(xlim[0], xlim[1], nx)
         y_list = numpy.linspace(ylim[0], ylim[1], ny)
         z_list = numpy.linspace(zlim[0], zlim[1], nz)
-        totalvolume = (xlim[1] - xlim[0]) * (ylim[1] - ylim[0]) * (zlim[1] - zlim[0])
+        totalvolume = abs(xlim[1] - xlim[0]) * abs(ylim[1] - ylim[0]) * abs(zlim[1] - zlim[0])
         vol = totalvolume / numberparticles
         if vol < 0:
             raise DomainError("Paritcles cannot have 0 volume")
@@ -964,7 +978,7 @@ def create_2D_domain(cls, xlim, ylim, nx, ny, type_id=1, mass=1.0,
         # Vertices
         x_list = numpy.linspace(xlim[0], xlim[1], nx)
         y_list = numpy.linspace(ylim[0], ylim[1], ny)
-        totalvolume = (xlim[1] - xlim[0]) * (ylim[1] - ylim[0])
+        totalvolume = abs(xlim[1] - xlim[0]) * abs(ylim[1] - ylim[0])
         vol = totalvolume / numberparticles
         if vol < 0:
             raise DomainError("Paritcles cannot have 0 volume")

spatialpy/core/domain.py

@@ -623,6 +623,8 @@ def validate(self, coverage="build", reactants=None, products=None, propensity_f
                         raise ReactionError("Type ids in restrict_to must be of type int or str.")
                     if type_id == "":
                         raise ReactionError("Type ids in restrict_to can't be an empty string.")
+                    if isinstance(type_id, str) and "UnAssigned" in type_id:
+                        raise ReactionError("'UnAssigned' is not a valid type_id.")
                     if coverage in ("all", "initialized"):
                         if not isinstance(type_id, str):
                             raise ReactionError("Type ids in restrict_to must be of type str.")

spatialpy/core/reaction.py

@@ -312,9 +312,9 @@ def __get_param_defs(self):
     def __get_particle_inits(self, num_chem_species):
         init_particles = ""
         if self.model.domain.type_id is None:
-            self.model.domain.type_id = ["type 1"] * self.model.domain.get_num_voxels()
+            self.model.domain.type_id = ["type_1"] * self.model.domain.get_num_voxels()
         for i, type_id in enumerate(self.model.domain.type_id):
-            if type_id is None:
+            if "UnAssigned" in type_id:
                 errmsg = "Not all particles have been defined in a type. Mass and other properties must be defined"
                 raise SimulationError(errmsg)
             x = self.model.domain.coordinates()[i, 0]

spatialpy/solvers/solver.py

@@ -173,12 +173,13 @@ def __build_element(stoich_species):
 def __get_particles(domain):
     s_particles = []
     for i, point in enumerate(domain.vertices):
+        type_id = domain.typeNdxMapping[domain.type_id[i]]
         s_particle = {"fixed":bool(domain.fixed[i]),
                       "mass":domain.mass[i],
                       "nu":domain.nu[i],
                       "particle_id":i,
                       "point":list(point),
-                      "type":int(domain.type_id[i]),
+                      "type":type_id,
                       "volume":domain.vol[i]}
 
         s_particles.append(s_particle)