diff --git a/pyproject.toml b/pyproject.toml index 470af73..4ed603d 100644 --- a/pyproject.toml +++ b/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 diff --git a/src/lineagetree/__init__.py b/src/lineagetree/__init__.py index 8c53c0e..b72ddcd 100644 --- a/src/lineagetree/__init__.py +++ b/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, diff --git a/src/lineagetree/measure/spatial.py b/src/lineagetree/measure/spatial.py index c74b55c..c7f4d8b 100644 --- a/src/lineagetree/measure/spatial.py +++ b/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 diff --git a/tests/test_lineageTree.py b/tests/test_lineageTree.py index 50bec7f..b8df55f 100644 --- a/tests/test_lineageTree.py +++ b/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():