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('_')] diff --git a/spatialpy/core/domain.py b/spatialpy/core/domain.py index f7afc481..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): """ @@ -474,7 +476,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 +485,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/lattice.py b/spatialpy/core/lattice.py index 2e7ca70a..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 @@ -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: 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()