HOTFIX: Lattice Numerical Error

spatialpy/core/__init__.py

@@ -18,33 +18,33 @@
 from spatialpy.__version__ import __version__
 from .boundarycondition import BoundaryCondition
 from .cleanup import (
-	cleanup_tempfiles,
-	cleanup_core_files,
-	cleanup_build_files,
-	cleanup_result_files
+    cleanup_tempfiles,
+    cleanup_core_files,
+    cleanup_build_files,
+    cleanup_result_files
 )
 from .datafunction import DataFunction
 from .domain import Domain
 from .geometry import (
-	CombinatoryGeometry,
-	Geometry,
-	GeometryAll,
-	GeometryExterior,
-	GeometryInterior
+    CombinatoryGeometry,
+    Geometry,
+    GeometryAll,
+    GeometryExterior,
+    GeometryInterior
 )
 from .initialcondition import (
-	InitialCondition,
-	PlaceInitialCondition,
-	UniformInitialCondition,
-	ScatterInitialCondition
+    InitialCondition,
+    PlaceInitialCondition,
+    UniformInitialCondition,
+    ScatterInitialCondition
 )
 from .lattice import (
-	CartesianLattice,
-	SphericalLattice,
-	CylindricalLattice,
-	XMLMeshLattice,
-	MeshIOLattice,
-	StochSSLattice
+    CartesianLattice,
+    SphericalLattice,
+    CylindricalLattice,
+    XMLMeshLattice,
+    MeshIOLattice,
+    StochSSLattice
 )
 from .model import Model, export_StochSS
 from .parameter import Parameter
@@ -54,21 +54,27 @@
 from .species import Species
 from .timespan import TimeSpan
 from .transformation import (
-	Transformation,
-	TranslationTransformation,
-	RotationTransformation,
-	ReflectionTransformation,
-	ScalingTransformation
+    Transformation,
+    TranslationTransformation,
+    RotationTransformation,
+    ReflectionTransformation,
+    ScalingTransformation
 )
 from .visualization import Visualization
 from .vtkreader import *
 
-_formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s')
-_handler = logging.StreamHandler()
-_handler.setFormatter(_formatter)
 version = __version__
-log = logging.getLogger()
+
+log = logging.getLogger("SpatialPy")
 log.setLevel(logging.WARNING)
-log.addHandler(_handler)
+log.propagate = False
+
+if not log.handlers:
+    _formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s')
+
+    _handler = logging.StreamHandler()
+    _handler.setFormatter(_formatter)
+
+    log.addHandler(_handler)
 
 __all__ = [s for s in dir() if not s.startswith('_')]

spatialpy/core/domain.py

@@ -474,7 +474,7 @@ def closest_vertex(self, point):
                 min_vtx = i
         return min_vtx
 
-    def compile_prep(self):
+    def compile_prep(self, allow_all_types=False):
         """
         Generate the domains list of type ids and check for invalid type_ids and rho values
         in preperation of compiling the simulation files.
@@ -483,8 +483,9 @@ def compile_prep(self):
         """
         self.apply_actions()
 
-        if self.type_id.tolist().count("type_UnAssigned") > 0:
-            raise DomainError("Particles must be assigned a type_id.")
+        if not allow_all_types:
+            if self.type_id.tolist().count("type_UnAssigned") > 0:
+                raise DomainError("Particles must be assigned a type_id.")
         if numpy.count_nonzero(self.rho) < len(self.rho):
             raise DomainError("Rho must be a positive value.")
 

spatialpy/core/lattice.py

@@ -158,8 +158,9 @@ def __add_point(self, domain, geometry, transform, point, count, kwargs):
             count += 1
         return count
 
-    def __generate_z(self, domain, geometry, transform, x, y, count, kwargs):
+    def __generate_z(self, domain, geometry, transform, x, y, count, z_digits, kwargs):
         for z in numpy.arange(self.zmin, self.zmax + self.deltaz, self.deltaz):
+            z = round(z, z_digits)
             count = self.__add_point(domain, geometry, transform, [x, y, z], count, kwargs)
         return count
 
@@ -190,14 +191,31 @@ def apply(self, domain, geometry, transform=None, **kwargs):
         if transform is not None and not callable(transform):
             raise LatticeError("transform must be a function.")
 
+        x_digits = max(
+            len(str(self.deltax).split(".")[1]) if "." in str(self.deltax) else 0,
+            len(str(self.xmin).split(".")[1]) if "." in str(self.xmin) else 0,
+            len(str(self.xmax).split(".")[1]) if "." in str(self.xmax) else 0
+        )
+        y_digits = max(
+            len(str(self.deltay).split(".")[1]) if "." in str(self.deltay) else 0,
+            len(str(self.ymin).split(".")[1]) if "." in str(self.ymin) else 0,
+            len(str(self.ymax).split(".")[1]) if "." in str(self.ymax) else 0
+        )
+        z_digits = max(
+            len(str(self.deltaz).split(".")[1]) if "." in str(self.deltaz) else 0,
+            len(str(self.zmin).split(".")[1]) if "." in str(self.zmin) else 0,
+            len(str(self.zmax).split(".")[1]) if "." in str(self.zmax) else 0
+        )
         count = 0
         for x in numpy.arange(self.xmin, self.xmax + self.deltax, self.deltax):
+            x = round(x, x_digits)
             for y in numpy.arange(self.ymin, self.ymax + self.deltay, self.deltay):
+                y = round(y, y_digits)
                 if self.deltaz == 0:
                     z = self.center[2]
                     count = self.__add_point(domain, geometry, transform, [x, y, z], count, kwargs)
                 else:
-                    count = self.__generate_z(domain, geometry, transform, x, y, count, kwargs)
+                    count = self.__generate_z(domain, geometry, transform, x, y, count, z_digits, kwargs)
         self._update_limits(domain)
         if 'vol' not in kwargs:
             offset = len(domain.vertices) - count
@@ -308,38 +326,48 @@ def apply(self, domain, geometry, transform=None, **kwargs):
         if transform is not None and not callable(transform):
             raise LatticeError("transform must be a function.")
 
+        digits = max(
+            len(str(self.deltar).split(".")[1]) if "." in str(self.deltar) else 0,
+            len(str(self.radius).split(".")[1]) if "." in str(self.radius) else 0
+        )
         count = 0
         radius = self.radius
         while radius > 0:
             # Calculate the approximate number of particle with the radius
             approx_rc = int(round((4 * radius ** 2) / ((self.deltas / 2) ** 2)))
 
-            # Set constants for the radius
-            p_area = 4 * numpy.pi * radius ** 2 / approx_rc
-            d_a = numpy.sqrt(p_area)
-            m_phi = int(round(numpy.pi * radius / d_a))
-            d_phi = numpy.pi / m_phi
-            d_theta = p_area / d_phi
-
-            for mphi in range(m_phi):
-                phi = numpy.pi * (mphi + 0.5) / m_phi
-                m_theta = int(round(2 * numpy.pi * numpy.sin(phi) / d_phi))
-
-                for mtheta in range(m_theta):
-                    theta = 2 * numpy.pi * mtheta / m_theta
-                    x = radius * numpy.cos(theta) * numpy.sin(phi)
-                    y = radius * numpy.sin(theta) * numpy.sin(phi)
-                    z = radius * numpy.cos(phi)
-                    if geometry.inside((x, y, z), False):
-                        if transform is None:
-                            point = [x, y, z]
-                        else:
-                            point = transform([x, y, z])
-                        if not isinstance(point, numpy.ndarray):
-                            point = numpy.array(point)
-                        domain.add_point(point + self.center, **kwargs)
-                        count += 1
-            radius -= self.deltar
+            if approx_rc == 0:
+                from spatialpy.core import log # pylint: disable=import-outside-toplevel
+                msg = f"Approximation of particles for the layer at radius {radius} is 0. "
+                msg += "Consider increasing the radius or increasing the radial spacing (deltas)"
+                log.warning(msg)
+            else:
+                # Set constants for the radius
+                p_area = 4 * numpy.pi * radius ** 2 / approx_rc
+                d_a = numpy.sqrt(p_area)
+                m_phi = int(round(numpy.pi * radius / d_a))
+                d_phi = numpy.pi / m_phi
+                d_theta = p_area / d_phi
+
+                for mphi in range(m_phi):
+                    phi = numpy.pi * (mphi + 0.5) / m_phi
+                    m_theta = int(round(2 * numpy.pi * numpy.sin(phi) / d_phi))
+
+                    for mtheta in range(m_theta):
+                        theta = 2 * numpy.pi * mtheta / m_theta
+                        x = radius * numpy.cos(theta) * numpy.sin(phi)
+                        y = radius * numpy.sin(theta) * numpy.sin(phi)
+                        z = radius * numpy.cos(phi)
+                        if geometry.inside((x, y, z), False):
+                            if transform is None:
+                                point = [x, y, z]
+                            else:
+                                point = transform([x, y, z])
+                            if not isinstance(point, numpy.ndarray):
+                                point = numpy.array(point)
+                            domain.add_point(point + self.center, **kwargs)
+                            count += 1
+            radius = round(radius - self.deltar, digits)
         if radius == 0 and geometry.inside((0, 0, 0), False):
             point = [0, 0, 0] if transform is None else transform([0, 0, 0])
             if not isinstance(point, numpy.ndarray):
@@ -439,6 +467,10 @@ def apply(self, domain, geometry, transform=None, **kwargs):
         if transform is not None and not callable(transform):
             raise LatticeError("transform must be a function.")
 
+        digits = max(
+            len(str(self.deltar).split(".")[1]) if "." in str(self.deltar) else 0,
+            len(str(self.radius).split(".")[1]) if "." in str(self.radius) else 0
+        )
         count = 0
         h_len = self.length / 2
         xmin = -h_len
@@ -448,28 +480,34 @@ def apply(self, domain, geometry, transform=None, **kwargs):
             # Calculate the approximate number of particle with the radius
             approx_rc = int(round((2 * radius * self.length) / ((self.deltas / 2) ** 2)))
 
-            p_area = 2 * numpy.pi * radius * self.length / approx_rc
-            d_a = numpy.sqrt(p_area)
-            m_theta = int(round(2 * numpy.pi * radius / d_a))
-            d_theta = 2 * numpy.pi / m_theta
-
-            x = xmin
-            while x <= xmax:
-                for mtheta in range(m_theta):
-                    theta = 2 * numpy.pi * (mtheta + 0.5) / m_theta
-                    y = radius * numpy.cos(theta)
-                    z = radius * numpy.sin(theta)
-                    if geometry.inside((x, y, z), False):
-                        if transform is None:
-                            point = [x, y, z]
-                        else:
-                            point = transform([x, y, z])
-                        if not isinstance(point, numpy.ndarray):
-                            point = numpy.array(point)
-                        domain.add_point(point + self.center, **kwargs)
-                        count += 1
-                x += self.deltas
-            radius -= self.deltar
+            if approx_rc == 0:
+                from spatialpy.core import log # pylint: disable=import-outside-toplevel
+                msg = f"Approximation of particles for the layer at radius {radius} is 0. "
+                msg += "Consider increasing the radius or increasing the radial spacing (deltas)"
+                log.warning(msg)
+            else:
+                p_area = 2 * numpy.pi * radius * self.length / approx_rc
+                d_a = numpy.sqrt(p_area)
+                m_theta = int(round(2 * numpy.pi * radius / d_a))
+                d_theta = 2 * numpy.pi / m_theta
+
+                x = xmin
+                while x <= xmax:
+                    for mtheta in range(m_theta):
+                        theta = 2 * numpy.pi * (mtheta + 0.5) / m_theta
+                        y = radius * numpy.cos(theta)
+                        z = radius * numpy.sin(theta)
+                        if geometry.inside((x, y, z), False):
+                            if transform is None:
+                                point = [x, y, z]
+                            else:
+                                point = transform([x, y, z])
+                            if not isinstance(point, numpy.ndarray):
+                                point = numpy.array(point)
+                            domain.add_point(point + self.center, **kwargs)
+                            count += 1
+                    x += self.deltas
+            radius = round(radius - self.deltar, digits)
         if radius == 0:
             x = xmin
             while x <= xmax:

spatialpy/core/model.py

@@ -447,7 +447,7 @@ def get_element(self, name):
             return self.get_data_function(name)
         raise ModelError(f"{self.name} does not contain an element named {name}.")
 
-    def add_domain(self, domain):
+    def add_domain(self, domain, allow_all_types=False):
         """
         Add a spatial domain to the model
 
@@ -461,7 +461,7 @@ def add_domain(self, domain):
                 "Unexpected parameter for add_domain. Parameter must be of type SpatialPy.Domain."
             )
 
-        domain.compile_prep()
+        domain.compile_prep(allow_all_types=allow_all_types)
         self.domain = domain
         return domain
 
@@ -766,7 +766,7 @@ def delete_all_reactions(self):
         self.listOfReactions.clear()
         self._listOfReactions.clear()
 
-    def get_reaction(self, rname):
+    def get_reaction(self, name):
         """
         Returns a reaction object by name.
 
@@ -979,7 +979,7 @@ def timespan(self, time_span, timestep_size=None):
         else:
             self.tspan = TimeSpan(time_span, timestep_size=timestep_size)
 
-    def compile_prep(self):
+    def compile_prep(self, allow_all_types=False):
         """
         Make sure all paramters are evaluated to scalars, update the models diffusion restrictions,
         create the models expression utility, and generate the domain list of type ids in preperation
@@ -997,7 +997,7 @@ def compile_prep(self):
 
         if self.domain is None:
             raise ModelError("The model's domain is not set.  Use 'add_domain()'.")
-        self.domain.compile_prep()
+        self.domain.compile_prep(allow_all_types=allow_all_types)
 
         self.__update_diffusion_restrictions()
         self.__apply_initial_conditions()