Adjusted to latest changes in py-pde (#55)

* Bumped py-pde required version to 0.35
* Removed the dependence on `grid.distance_real` and
`grid.polar_coordinates_real`

docs/source/installation.rst

@@ -41,7 +41,7 @@ matplotlib   >=3.1.0   Visualizing results
 numpy        >=1.22    Array library used for storing data
 numba        >=0.48    Just-in-time compilation to accelerate numerics
 scipy        >=1.4     Miscellaneous scientific functions
-py-pde       >=0.34    Simulating partial differential equations
+py-pde       >=0.35    Simulating partial differential equations
 ===========  ========= =========
 
 These package can be installed via your operating system's package manager, e.g.

droplets/droplet_tracks.py

@@ -7,12 +7,12 @@
    DropletTrack
    DropletTrackList
 
-
 .. codeauthor:: David Zwicker <[email protected]>
 """
 
 from __future__ import annotations
 
+import functools
 import json
 import logging
 from typing import Callable, Literal
@@ -448,7 +448,7 @@ def plot_positions(
             segments = []
             for p1, p2 in zip(xy[:-1], xy[1:]):
                 dist_direct = np.hypot(p1[0] - p2[0], p1[1] - p2[1])
-                dist_real = grid.distance_real(p1, p2)
+                dist_real = grid.distance(p1, p2, coords="cartesian")
                 close = np.isclose(dist_direct, dist_real)
                 segments.append(close)
 
@@ -568,7 +568,7 @@ def match_tracks(
                     if grid is None:
                         metric: str | Callable = "euclidean"
                     else:
-                        metric = grid.distance_real
+                        metric = functools.partial(grid.distance, coords="cartesian")
                     points_prev = [track.last.position for track in tracks_alive]
                     points_now = [droplet.position for droplet in emulsion]
                     dists = distance.cdist(points_prev, points_now, metric=metric)

droplets/droplets.py

@@ -32,7 +32,7 @@
 import warnings
 from abc import ABCMeta, abstractmethod
 from pathlib import Path
-from typing import Any, Callable, List, Tuple, Type, TypeVar, Union
+from typing import Any, Callable, List, Literal, Tuple, Type, TypeVar, Union, overload
 
 import numpy as np
 from numba.extending import register_jitable
@@ -41,7 +41,7 @@
 from scipy import integrate
 
 from pde.fields import ScalarField
-from pde.grids.base import GridBase
+from pde.grids.base import DimensionError, GridBase
 from pde.tools.cuboid import Cuboid
 from pde.tools.misc import preserve_scalars
 from pde.tools.plotting import PlotReference, plot_on_axes
@@ -96,6 +96,77 @@ def iterate_in_pairs(it, fill=0):
             break
 
 
+@overload
+def polar_coordinates(
+    grid: GridBase,
+    *,
+    origin: np.ndarray | None = None,
+    ret_angle: Literal[False] = False,
+) -> np.ndarray:
+    ...
+
+
+@overload
+def polar_coordinates(
+    grid: GridBase, *, origin: np.ndarray | None = None, ret_angle: Literal[True]
+) -> tuple[np.ndarray, ...]:
+    ...
+
+
+def polar_coordinates(
+    grid: GridBase, *, origin: np.ndarray | None = None, ret_angle: bool = False
+) -> np.ndarray | tuple[np.ndarray, ...]:
+    """return polar coordinates associated with grid points
+
+    Args:
+        grid (:class:`~pde.grids.base.GridBase`):
+            The grid whose cell coordinates are used.
+        origin (:class:`~numpy.ndarray`, optional):
+            Cartesian coordinates of the origin at which polar coordinates are anchored.
+        ret_angle (bool):
+            Determines whether angles are returned alongside the distance. If `False`
+            only the distance to the origin is returned for each support point of the
+            grid. If `True`, the distance and angles are returned. For a 1d system
+            system, the angle is defined as the sign of the difference between the
+            point and the origin, so that angles can either be 1 or -1. For 2d
+            systems and 3d systems, polar coordinates and spherical coordinates are
+            used, respectively.
+
+    Returns:
+        :class:`~numpy.ndarray` or tuple of :class:`~numpy.ndarray`:
+            Coordinates values in polar coordinates
+    """
+    if origin is None:
+        origin = np.zeros(grid.dim)
+    else:
+        origin = np.asarray(origin, dtype=float)
+        if origin.shape != (grid.dim,):
+            raise DimensionError("Dimensions are not compatible")
+
+    # calculate the difference vector between all cells and the origin
+    origin_grid = grid.transform(origin, source="cartesian", target="grid")
+    diff = grid.difference_vector(grid.cell_coords, origin_grid)
+    dist: np.ndarray = np.linalg.norm(diff, axis=-1)  # get distance
+
+    # determine distance and optionally angles for these vectors
+    if not ret_angle:
+        return dist
+
+    elif grid.dim == 1:
+        return dist, np.sign(diff)[..., 0]
+
+    elif grid.dim == 2:
+        return dist, np.arctan2(diff[..., 0], diff[..., 1])
+
+    elif grid.dim == 3:
+        theta = np.arccos(diff[..., 2] / dist)
+        phi = np.arctan2(diff[..., 0], diff[..., 1])
+        return dist, theta, phi
+
+    else:
+        raise NotImplementedError(f"Cannot calculate angles for dimension {grid.dim}")
+
+
 class DropletBase:
     """represents a generic droplet
 
@@ -376,7 +447,7 @@ def overlaps(self, other: SphericalDroplet, grid: GridBase | None = None) -> boo
         if grid is None:
             distance = float(np.linalg.norm(self.position - other.position))
         else:
-            distance = grid.distance_real(self.position, other.position)
+            distance = grid.distance(self.position, other.position, coords="cartesian")
         return distance < self.radius + other.radius
 
     @classmethod
@@ -464,8 +535,8 @@ def _get_phase_field(self, grid: GridBase, dtype: DTypeLike = float) -> np.ndarr
             )
 
         # calculate distances from droplet center
-        dist = grid.polar_coordinates_real(self.position)
-        return (dist < self.radius).astype(dtype)  # type: ignore
+        dist = polar_coordinates(grid, origin=self.position, ret_angle=False)
+        return (dist < self.radius).astype(dtype)
 
     def get_phase_field(
         self,
@@ -692,7 +763,7 @@ def _get_phase_field(self, grid: GridBase, dtype: DTypeLike = float) -> np.ndarr
             interface_width = self.interface_width
 
         # calculate distances from droplet center
-        dist: np.ndarray = grid.polar_coordinates_real(self.position, ret_angle=False)  # type: ignore
+        dist = polar_coordinates(grid, origin=self.position, ret_angle=False)
 
         # make the image
         if interface_width == 0 or np.issubdtype(dtype, bool):
@@ -838,7 +909,7 @@ def _get_phase_field(self, grid: GridBase, dtype: DTypeLike = float) -> np.ndarr
             interface_width = self.interface_width
 
         # calculate grid distance from droplet center
-        dist, *angles = grid.polar_coordinates_real(self.position, ret_angle=True)
+        dist, *angles = polar_coordinates(grid, origin=self.position, ret_angle=True)
 
         # calculate interface distance from droplet center
         interface = self.interface_distance(*angles)

droplets/emulsions.py

@@ -502,7 +502,7 @@ def get_distance(p1, p2):
                 return np.linalg.norm(p1 - p2)
 
         else:
-            get_distance = grid.distance_real
+            get_distance = functools.partial(grid.distance, coords="cartesian")
 
         # calculate pairwise distance and return it in requested form
         num = len(self)

pyproject.toml

@@ -30,7 +30,7 @@ dependencies = [
     "numpy>=1.22.0",
     "scipy>=1.4.0",
     "sympy>=1.5.0",
-    "py-pde>=0.34",
+    "py-pde>=0.35",
 ]
 
 [project.optional-dependencies]

requirements.txt

@@ -2,4 +2,4 @@ matplotlib>=3.1.0
 numpy>=1.22.0
 numba>=0.48.0
 scipy>=1.4.0
-py-pde>=0.34
+py-pde>=0.35

scripts/run_tests.py

@@ -155,7 +155,14 @@ def run_unit_tests(
     args.append("tests")
 
     # actually run the test
-    return sp.run(args, env=env, cwd=PACKAGE_PATH).returncode
+    retcode = sp.run(args, env=env, cwd=PACKAGE_PATH).returncode
+
+    # delete intermediate coverage files, which are sometimes left behind
+    if coverage:
+        for p in Path("..").glob(".coverage*"):
+            p.unlink()
+
+    return retcode
 
 
 def main() -> int: