From b7875b3536f44c542c71f42fd01e4ae901ae099c Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Wed, 25 Jan 2023 10:28:00 -0500 Subject: [PATCH 1/9] fixed numerical error when decrementing radius. --- spatialpy/core/lattice.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/spatialpy/core/lattice.py b/spatialpy/core/lattice.py index 2e7ca70a..cc573aa3 100644 --- a/spatialpy/core/lattice.py +++ b/spatialpy/core/lattice.py @@ -308,6 +308,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 radius = self.radius while radius > 0: @@ -323,7 +327,7 @@ def apply(self, domain, geometry, transform=None, **kwargs): 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)) + m_theta = int(round(2 * numpy.pi * numpy.sin(phi) / d_theta)) for mtheta in range(m_theta): theta = 2 * numpy.pi * mtheta / m_theta @@ -339,7 +343,7 @@ def apply(self, domain, geometry, transform=None, **kwargs): point = numpy.array(point) domain.add_point(point + self.center, **kwargs) count += 1 - radius -= self.deltar + 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 +443,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 @@ -456,7 +464,7 @@ def apply(self, domain, geometry, transform=None, **kwargs): x = xmin while x <= xmax: for mtheta in range(m_theta): - theta = 2 * numpy.pi * (mtheta + 0.5) / m_theta + theta = 2 * numpy.pi * (mtheta + 0.5) / d_theta y = radius * numpy.cos(theta) z = radius * numpy.sin(theta) if geometry.inside((x, y, z), False): @@ -469,7 +477,7 @@ def apply(self, domain, geometry, transform=None, **kwargs): domain.add_point(point + self.center, **kwargs) count += 1 x += self.deltas - radius -= self.deltar + radius = round(radius - self.deltar, digits) if radius == 0: x = xmin while x <= xmax: From bbdc022a5b0c93cf845458726bad41d766b5f3f8 Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Wed, 25 Jan 2023 10:59:28 -0500 Subject: [PATCH 2/9] fixed numerical error when decrementing radius. --- spatialpy/core/lattice.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/spatialpy/core/lattice.py b/spatialpy/core/lattice.py index cc573aa3..6d1b43d9 100644 --- a/spatialpy/core/lattice.py +++ b/spatialpy/core/lattice.py @@ -327,7 +327,7 @@ def apply(self, domain, geometry, transform=None, **kwargs): 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_theta)) + m_theta = int(round(2 * numpy.pi * numpy.sin(phi) / m_theta)) for mtheta in range(m_theta): theta = 2 * numpy.pi * mtheta / m_theta @@ -464,7 +464,7 @@ def apply(self, domain, geometry, transform=None, **kwargs): x = xmin while x <= xmax: for mtheta in range(m_theta): - theta = 2 * numpy.pi * (mtheta + 0.5) / d_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): From f3af2a4042b2210f1eecde77b7ec06ad96663a97 Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Wed, 25 Jan 2023 11:30:31 -0500 Subject: [PATCH 3/9] fixed numerical error when decrementing radius. --- spatialpy/core/lattice.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/spatialpy/core/lattice.py b/spatialpy/core/lattice.py index 6d1b43d9..dd9e855c 100644 --- a/spatialpy/core/lattice.py +++ b/spatialpy/core/lattice.py @@ -327,7 +327,7 @@ def apply(self, domain, geometry, transform=None, **kwargs): for mphi in range(m_phi): phi = numpy.pi * (mphi + 0.5) / m_phi - m_theta = int(round(2 * numpy.pi * numpy.sin(phi) / m_theta)) + 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 From 60a87c8a6f030d34f244df7cefe7dd13bf36c188 Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Wed, 25 Jan 2023 13:22:18 -0500 Subject: [PATCH 4/9] Added warning for approx_rc == 0. --- spatialpy/core/lattice.py | 104 +++++++++++++++++++++----------------- 1 file changed, 58 insertions(+), 46 deletions(-) diff --git a/spatialpy/core/lattice.py b/spatialpy/core/lattice.py index dd9e855c..9724ad62 100644 --- a/spatialpy/core/lattice.py +++ b/spatialpy/core/lattice.py @@ -318,31 +318,37 @@ def apply(self, domain, geometry, transform=None, **kwargs): # 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 + 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]) @@ -456,27 +462,33 @@ 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 + 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 From 1e780930191b461de461267d2f7e8bf76ec12f23 Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Mon, 30 Jan 2023 09:36:42 -0500 Subject: [PATCH 5/9] Changed the logger to a named logger. --- spatialpy/core/__init__.py | 64 +++++++++++++++++++++----------------- 1 file changed, 35 insertions(+), 29 deletions(-) diff --git a/spatialpy/core/__init__.py b/spatialpy/core/__init__.py index af272085..af5e9d71 100644 --- a/spatialpy/core/__init__.py +++ b/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('_')] From 8dcd124b25ab6af03db7848740d19405f71445a7 Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Mon, 30 Jan 2023 10:20:28 -0500 Subject: [PATCH 6/9] Fixed warning messages for approx_rc value of 0. --- spatialpy/core/lattice.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/spatialpy/core/lattice.py b/spatialpy/core/lattice.py index 9724ad62..2164adb2 100644 --- a/spatialpy/core/lattice.py +++ b/spatialpy/core/lattice.py @@ -320,7 +320,7 @@ def apply(self, domain, geometry, transform=None, **kwargs): 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 = 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: @@ -464,7 +464,7 @@ def apply(self, domain, geometry, transform=None, **kwargs): 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 = 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: From c9673b401830215327166679a7dcd3928e94f222 Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Mon, 30 Jan 2023 14:37:23 -0500 Subject: [PATCH 7/9] Added flag to allow all types when executing compile prep. --- spatialpy/core/domain.py | 7 ++++--- spatialpy/core/model.py | 10 +++++----- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/spatialpy/core/domain.py b/spatialpy/core/domain.py index f7afc481..754a99b3 100644 --- a/spatialpy/core/domain.py +++ b/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.") diff --git a/spatialpy/core/model.py b/spatialpy/core/model.py index 25353e31..f5fe3687 100644 --- a/spatialpy/core/model.py +++ b/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() From 72ec726590dacc049fd9a06e7a6183778d97f3f0 Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Mon, 30 Jan 2023 15:14:44 -0500 Subject: [PATCH 8/9] Added handler for math errors for cartesian lattices. --- spatialpy/core/lattice.py | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/spatialpy/core/lattice.py b/spatialpy/core/lattice.py index 2164adb2..b9c5f5f6 100644 --- a/spatialpy/core/lattice.py +++ b/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 From 6d56050f7a2e673bbf4163ddfd26eefc3ccf66cd Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Tue, 31 Jan 2023 15:27:36 -0500 Subject: [PATCH 9/9] Fixed issues with removing multiple particles from the domain. --- spatialpy/core/domain.py | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/spatialpy/core/domain.py b/spatialpy/core/domain.py index 754a99b3..8bbf22ec 100644 --- a/spatialpy/core/domain.py +++ b/spatialpy/core/domain.py @@ -388,17 +388,19 @@ def apply_remove_action(self, action): # remove the particles that fall within the defined region on_boundary = self.find_boundary_points(update=True) - for v_ndx, vertex in enumerate(self.vertices): - if action['geometry'].inside(vertex, on_boundary[v_ndx]): - self.vertices = numpy.delete(self.vertices, v_ndx, 0) - self.type_id = numpy.delete(self.type_id, v_ndx, 0) - self.vol = numpy.delete(self.vol, v_ndx) - self.mass = numpy.delete(self.mass, v_ndx) - self.rho = numpy.delete(self.rho, v_ndx) - self.nu = numpy.delete(self.nu, v_ndx) - self.c = numpy.delete(self.c, v_ndx) - self.fixed = numpy.delete(self.fixed, v_ndx) - self.on_boundary = numpy.delete(self.on_boundary, v_ndx) + geometry = action['geometry'] + indices = [v_ndx for v_ndx, vertex in enumerate(self.vertices) if geometry.inside(vertex, on_boundary[v_ndx])] + indices.reverse() + for v_ndx in indices: + self.vertices = numpy.delete(self.vertices, v_ndx, 0) + self.type_id = numpy.delete(self.type_id, v_ndx, 0) + self.vol = numpy.delete(self.vol, v_ndx) + self.mass = numpy.delete(self.mass, v_ndx) + self.rho = numpy.delete(self.rho, v_ndx) + self.nu = numpy.delete(self.nu, v_ndx) + self.c = numpy.delete(self.c, v_ndx) + self.fixed = numpy.delete(self.fixed, v_ndx) + self.on_boundary = numpy.delete(self.on_boundary, v_ndx) def apply_set_action(self, action): """