fix for all timepoints

pyproject.toml

@@ -16,7 +16,7 @@ maintainers = [
 ]
 name = "lineagetree"
 description = "Structure for Lineage Trees"
-version = "3.0.4"
+version = "3.1.0"
 license = "MIT"
 license-files = [ "LICENSE" ]
 readme = {file = "README.md", content-type = "text/markdown"}
@@ -76,7 +76,7 @@ profile = "black"
 line_length = 79
 
 [tool.bumpver]
-current_version = "3.0.4"
+current_version = "3.1.0"
 version_pattern = "MAJOR.MINOR.PATCH[-TAG]"
 commit_message = "bump version {old_version} -> {new_version}"
 commit = true

src/lineagetree/__init__.py

@@ -1,4 +1,4 @@
-__version__ = "3.0.4"
+__version__ = "3.1.0"
 from .lineage_tree import LineageTree
 from ._io._loaders import (
     read_from_ASTEC,

src/lineagetree/measure/spatial.py

@@ -1,7 +1,8 @@
 from __future__ import annotations
+from warnings import warn, catch_warnings, simplefilter
 
 from itertools import combinations
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, Iterable
 
 import numpy as np
 from scipy.spatial import Delaunay, KDTree
@@ -49,7 +50,9 @@ def get_idx3d(lT: LineageTree, t: int) -> tuple[KDTree, np.ndarray]:
     return idx3d, np.array(to_check_lT)
 
 
-def get_gabriel_graph(lT: LineageTree, t: int) -> dict[int, set[int]]:
+def get_gabriel_graph(
+    lT: LineageTree, time: int | Iterable[int] | None = None
+) -> dict[int, set[int]]:
     """Build the Gabriel graph of the given graph for time point `t`.
     The Garbiel graph is then stored in `lT.Gabriel_graph` and returned.
 
@@ -59,8 +62,9 @@ def get_gabriel_graph(lT: LineageTree, t: int) -> dict[int, set[int]]:
     ----------
     lT : LineageTree
         The LineageTree instance.
-    t : int
-        time
+    time : int or Iterable of int, optional
+        time or iterable of times.
+        If not given the gabriel graph will be calculated for all timepoints.
 
     Returns
     -------
@@ -70,43 +74,56 @@ def get_gabriel_graph(lT: LineageTree, t: int) -> dict[int, set[int]]:
     if not hasattr(lT, "Gabriel_graph"):
         lT.Gabriel_graph = {}
 
-    if lT.time_nodes[t] - lT.Gabriel_graph.keys():
-        _, nodes = lT.get_idx3d(t)
-
-        data_corres = {}
-        data = []
-        for i, C in enumerate(nodes):
-            data.append(lT.pos[C])
-            data_corres[i] = C
-
-        tmp = Delaunay(data)
-
-        delaunay_graph = {}
-
-        for N in tmp.simplices:
-            for e1, e2 in combinations(np.sort(N), 2):
-                delaunay_graph.setdefault(e1, set()).add(e2)
-                delaunay_graph.setdefault(e2, set()).add(e1)
-
-        Gabriel_graph = {}
-
-        for e1, neighbs in delaunay_graph.items():
-            for ni in neighbs:
-                if not any(
-                    np.linalg.norm((data[ni] + data[e1]) / 2 - data[i])
-                    < np.linalg.norm(data[ni] - data[e1]) / 2
-                    for i in delaunay_graph[e1].intersection(
-                        delaunay_graph[ni]
-                    )
-                ):
-                    Gabriel_graph.setdefault(data_corres[e1], set()).add(
-                        data_corres[ni]
-                    )
-                    Gabriel_graph.setdefault(data_corres[ni], set()).add(
-                        data_corres[e1]
-                    )
-
-        lT.Gabriel_graph.update(Gabriel_graph)
+    if time is None:
+        time = lT.time_nodes.keys()
+    elif not isinstance(time, Iterable):
+        time = [time]
+
+    for t in time:
+        if lT.time_nodes[t] - lT.Gabriel_graph.keys():
+            nodes = np.fromiter(list(lT.time_nodes[t]), dtype=int)
+
+            data_corres = {}
+            data = []
+            for i, C in enumerate(nodes):
+                data.append(lT.pos[C])
+                data_corres[i] = C
+
+            delaunay_graph = {}
+
+            # The delaunay triangulation is only usefult to compute
+            # when the number of points is higher than the spatial dimension + 1
+            if len(data[0]) + 1 < len(data):
+                tmp = Delaunay(data)
+                for N in tmp.simplices:
+                    for e1, e2 in combinations(np.sort(N), 2):
+                        delaunay_graph.setdefault(e1, set()).add(e2)
+                        delaunay_graph.setdefault(e2, set()).add(e1)
+            # When there are fewer nodes than the number of dimensions + 2
+            # The Delaunay is the complete graph
+            else:
+                for e1, e2 in combinations(data_corres, 2):
+                    delaunay_graph.setdefault(e1, set()).add(e2)
+                    delaunay_graph.setdefault(e2, set()).add(e1)
+
+            Gabriel_graph = {}
+
+            for e1, neighbs in delaunay_graph.items():
+                for ni in neighbs:
+                    if not any(
+                        np.linalg.norm((data[ni] + data[e1]) / 2 - data[i])
+                        < np.linalg.norm(data[ni] - data[e1]) / 2
+                        for i in delaunay_graph[e1].intersection(
+                            delaunay_graph[ni]
+                        )
+                    ):
+                        Gabriel_graph.setdefault(data_corres[e1], set()).add(
+                            data_corres[ni]
+                        )
+                        Gabriel_graph.setdefault(data_corres[ni], set()).add(
+                            data_corres[e1]
+                        )
+            lT.Gabriel_graph.update(Gabriel_graph)
 
     return lT.Gabriel_graph
 

tests/test_lineageTree.py

@@ -585,6 +585,11 @@ def test_idx3d():
 def test_gabriel_graph():
     gg = lT1.get_gabriel_graph(0)
     assert gg[173618] == {110832, 168322}
+    gg_all = lT1.get_gabriel_graph()
+    gg_all_2 = {}
+    for t in lT1.time_nodes:
+        gg_all_2.update(lT1.get_gabriel_graph(t))
+    assert gg_all == gg_all_2
 
 
 def test_get_chain_of_node():