From 78699ec220f1cd1fa19b1943b48bb1f133ac045c Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 17 Jun 2022 15:10:34 +0300 Subject: [PATCH] add recursive clustering and skeletonization --- doc/conf.py | 4 + doc/linalg.rst | 22 ++- pytential/linalg/__init__.py | 16 -- pytential/linalg/cluster.py | 277 ++++++++++++++++++++++++++++ pytential/linalg/proxy.py | 72 +------- pytential/linalg/skeletonization.py | 176 +++++++++++++----- pytential/linalg/utils.py | 8 +- pytential/symbolic/matrix.py | 5 +- test/extra_matrix_data.py | 10 +- test/test_linalg_cluster.py | 131 +++++++++++++ test/test_linalg_proxy.py | 10 +- test/test_linalg_skeletonization.py | 57 ++++-- test/test_matrix.py | 6 +- 13 files changed, 619 insertions(+), 175 deletions(-) create mode 100644 pytential/linalg/cluster.py create mode 100644 test/test_linalg_cluster.py diff --git a/doc/conf.py b/doc/conf.py index 636bb774b..38a0ead35 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -22,6 +22,10 @@ "DOFDescriptorLike": "pytential.symbolic.dof_desc.DOFDescriptorLike", } +nitpick_ignore_regex = [ + ["py:class", r"_ProxyNeighborEvaluationResult"], + ] + intersphinx_mapping = { "https://docs.python.org/3/": None, "https://numpy.org/doc/stable/": None, diff --git a/doc/linalg.rst b/doc/linalg.rst index b75f7dc5c..0a9b6ccb0 100644 --- a/doc/linalg.rst +++ b/doc/linalg.rst @@ -8,24 +8,38 @@ scheme is used: component of the Stokeslet. * ``cluster`` refers to a piece of a ``block`` as used by the recursive proxy-based skeletonization of the direct solver algorithms. Clusters - are represented by a :class:`~pytential.linalg.TargetAndSourceClusterList`. + are represented by a :class:`~pytential.linalg.utils.TargetAndSourceClusterList`. .. _direct_solver: Hierarchical Direct Solver -------------------------- +.. note:: + + High-level API for direct solvers is in progress. + +Low-level Functionality +----------------------- + .. warning:: All the classes and routines in this module are experimental and the API can change at any point. +.. automodule:: pytential.linalg.skeletonization +.. automodule:: pytential.linalg.cluster .. automodule:: pytential.linalg.proxy -.. automodule:: pytential.linalg.utils -Internal Functionality ----------------------- +Internal Functionality and Utilities +------------------------------------ + +.. warning:: + All the classes and routines in this module are experimental and the + API can change at any point. + +.. automodule:: pytential.linalg.utils .. automodule:: pytential.linalg.direct_solver_symbolic .. vim: sw=4:tw=75 diff --git a/pytential/linalg/__init__.py b/pytential/linalg/__init__.py index bfbd4d82e..eed539cc8 100644 --- a/pytential/linalg/__init__.py +++ b/pytential/linalg/__init__.py @@ -25,25 +25,9 @@ make_index_list, make_index_cluster_cartesian_product, interp_decomp, ) -from pytential.linalg.proxy import ( - ProxyClusterGeometryData, ProxyPointTarget, ProxyPointSource, - ProxyGeneratorBase, ProxyGenerator, QBXProxyGenerator, - partition_by_nodes, gather_cluster_neighbor_points, - ) -from pytential.linalg.skeletonization import ( - SkeletonizationWrangler, make_skeletonization_wrangler, - SkeletonizationResult, skeletonize_by_proxy, - ) __all__ = ( "IndexList", "TargetAndSourceClusterList", "make_index_list", "make_index_cluster_cartesian_product", "interp_decomp", - - "ProxyClusterGeometryData", "ProxyPointTarget", "ProxyPointSource", - "ProxyGeneratorBase", "ProxyGenerator", "QBXProxyGenerator", - "partition_by_nodes", "gather_cluster_neighbor_points", - - "SkeletonizationWrangler", "make_skeletonization_wrangler", - "SkeletonizationResult", "skeletonize_by_proxy", ) diff --git a/pytential/linalg/cluster.py b/pytential/linalg/cluster.py new file mode 100644 index 000000000..21d9f3eda --- /dev/null +++ b/pytential/linalg/cluster.py @@ -0,0 +1,277 @@ +__copyright__ = "Copyright (C) 2022 Alexandru Fikl" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +from dataclasses import dataclass, replace +from functools import singledispatch +from typing import Optional + +import numpy as np + +from pytools import T, memoize_method + +from arraycontext import PyOpenCLArrayContext +from boxtree.tree import Tree +from pytential import sym, GeometryCollection +from pytential.linalg.utils import IndexList, TargetAndSourceClusterList + +__doc__ = """ +Clustering +~~~~~~~~~~ + +.. autoclass:: ClusterTreeLevel + +.. autofunction:: cluster +.. autofunction:: partition_by_nodes +""" + +# FIXME: this is just an arbitrary value +_DEFAULT_MAX_PARTICLES_IN_BOX = 32 + + +# {{{ cluster tree + +@dataclass(frozen=True) +class ClusterTreeLevel: + """ + .. attribute:: level + + Current level that is represented. + + .. attribute:: nlevels + + Total number of levels in the tree. + + .. attribute:: nclusters + + Number of clusters on the current level (same as number of boxes + in :attr:`partition_box_ids`). + + .. attribute:: partition_box_ids + + Box IDs on the current level. + + .. attribute:: partition_parent_ids + + Parent box IDs for :attr:`partition_box_ids`. + + .. attribute:: partition_parent_map + + An object :class:`~numpy.ndarray`, where each entry maps a parent + to its children indices in :attr:`partition_box_ids`. This can be used to + :func:`cluster` all indices in ``partition_parent_map[i]`` into their + parent. + + .. automethod:: parent + """ + + # level info + level: int + partition_box_ids: np.ndarray + + # tree info + nlevels: int + box_parent_ids: np.ndarray + + # NOTE: only here to allow easier debugging + testing + _tree: Optional[Tree] + + @property + def nclusters(self): + return self.partition_box_ids.size + + @property + def partition_parent_ids(self): + return self.box_parent_ids[self.partition_box_ids] + + @property + @memoize_method + def partition_parent_map(self): + # NOTE: np.unique returns a sorted array + unique_parent_ids = np.unique(self.partition_parent_ids) + # find the index of each parent id + unique_parent_index = np.searchsorted( + unique_parent_ids, self.partition_parent_ids + ) + + unique_parent_map = np.empty(unique_parent_ids.size, dtype=object) + for i in range(unique_parent_ids.size): + unique_parent_map[i], = np.nonzero(unique_parent_index == i) + + return unique_parent_map + + def parent(self) -> "ClusterTreeLevel": + """ + :returns: a new :class:`ClusterTreeLevel` that represents the parent of + the current one, with appropriately updated :attr:`partition_box_ids`, + etc. + """ + + if self.nclusters == 1: + assert self.level == 0 + return self + + return replace(self, + level=self.level - 1, + partition_box_ids=np.unique(self.partition_parent_ids)) + + +@singledispatch +def cluster(obj: T, ctree: ClusterTreeLevel) -> T: + """Merge together elements of *obj* into their parent object, as described + by *ctree*. + """ + raise NotImplementedError(type(obj).__name__) + + +@cluster.register(IndexList) +def _cluster_index_list(obj: IndexList, ctree: ClusterTreeLevel) -> IndexList: + assert obj.nclusters == ctree.nclusters + + if ctree.nclusters == 1: + return obj + + sizes = [ + sum(obj.cluster_size(j) for j in ppm) + for ppm in ctree.partition_parent_map + ] + return replace(obj, starts=np.cumsum([0] + sizes)) + + +@cluster.register(TargetAndSourceClusterList) +def _cluster_target_and_source_cluster_list( + obj: TargetAndSourceClusterList, ctree: ClusterTreeLevel, + ) -> TargetAndSourceClusterList: + assert obj.nclusters == ctree.nclusters + + if ctree.nclusters == 1: + return obj + + return replace(obj, + targets=cluster(obj.targets, ctree), + sources=cluster(obj.sources, ctree)) + +# }}} + + +# {{{ cluster generation + +def _build_binary_ish_tree_from_indices(starts: np.ndarray) -> ClusterTreeLevel: + partition_box_ids = np.arange(starts.size - 1) + + box_ids = partition_box_ids + + box_parent_ids = [] + offset = box_ids.size + while box_ids.size > 1: + # NOTE: this is probably not the most efficient way to do it, but this + # code is mostly meant for debugging using a simple tree + clusters = np.array_split(box_ids, box_ids.size // 2) + parent_ids = offset + np.arange(len(clusters)) + box_parent_ids.append(np.repeat(parent_ids, [len(c) for c in clusters])) + + box_ids = parent_ids + offset += box_ids.size + + # NOTE: make the root point to itself + box_parent_ids.append(np.array([offset - 1])) + nlevels = len(box_parent_ids) + + return ClusterTreeLevel( + level=nlevels - 1, + partition_box_ids=partition_box_ids, + nlevels=nlevels, + box_parent_ids=np.concatenate(box_parent_ids), + _tree=None) + + +def partition_by_nodes( + actx: PyOpenCLArrayContext, places: GeometryCollection, *, + dofdesc: Optional[sym.DOFDescriptorLike] = None, + tree_kind: Optional[str] = "adaptive-level-restricted", + max_particles_in_box: Optional[int] = None) -> IndexList: + """Generate equally sized ranges of nodes. The partition is created at the + lowest level of granularity, i.e. nodes. This results in balanced ranges + of points, but will split elements across different ranges. + + :arg dofdesc: a :class:`~pytential.symbolic.dof_desc.DOFDescriptor` for + the geometry in *places* which should be partitioned. + :arg tree_kind: if not *None*, it is passed to :class:`boxtree.TreeBuilder`. + :arg max_particles_in_box: value used to control the number of points + in each partition (and thus the number of partitions). See the documentation + in :class:`boxtree.TreeBuilder`. + """ + if dofdesc is None: + dofdesc = places.auto_source + dofdesc = sym.as_dofdesc(dofdesc) + + if max_particles_in_box is None: + max_particles_in_box = _DEFAULT_MAX_PARTICLES_IN_BOX + + lpot_source = places.get_geometry(dofdesc.geometry) + discr = places.get_discretization(dofdesc.geometry, dofdesc.discr_stage) + + if tree_kind is not None: + from pytential.qbx.utils import tree_code_container + tcc = tree_code_container(lpot_source._setup_actx) + + from arraycontext import flatten + from meshmode.dof_array import DOFArray + tree, _ = tcc.build_tree()(actx.queue, + particles=flatten( + actx.thaw(discr.nodes()), actx, leaf_class=DOFArray + ), + max_particles_in_box=max_particles_in_box, + kind=tree_kind) + + from boxtree import box_flags_enum + tree = tree.get(actx.queue) + leaf_boxes, = (tree.box_flags & box_flags_enum.HAS_CHILDREN == 0).nonzero() + + indices = np.empty(len(leaf_boxes), dtype=object) + starts = None + + for i, ibox in enumerate(leaf_boxes): + box_start = tree.box_source_starts[ibox] + box_end = box_start + tree.box_source_counts_cumul[ibox] + indices[i] = tree.user_source_ids[box_start:box_end] + + ctree = ClusterTreeLevel( + level=tree.nlevels - 1, + nlevels=tree.nlevels, + box_parent_ids=tree.box_parent_ids, + partition_box_ids=leaf_boxes, + _tree=tree) + else: + if discr.ambient_dim != 2 and discr.dim == 1: + raise ValueError("only curves are supported for 'tree_kind=None'") + + nclusters = max(discr.ndofs // max_particles_in_box, 2) + indices = np.arange(0, discr.ndofs, dtype=np.int64) + starts = np.linspace(0, discr.ndofs, nclusters + 1, dtype=np.int64) + assert starts[-1] == discr.ndofs + + ctree = _build_binary_ish_tree_from_indices(starts) + + from pytential.linalg import make_index_list + return make_index_list(indices, starts=starts), ctree + +# }}} diff --git a/pytential/linalg/proxy.py b/pytential/linalg/proxy.py index 642879316..37c907728 100644 --- a/pytential/linalg/proxy.py +++ b/pytential/linalg/proxy.py @@ -45,8 +45,6 @@ Proxy Point Generation ~~~~~~~~~~~~~~~~~~~~~~ -.. currentmodule:: pytential.linalg - .. autoclass:: ProxyPointSource .. autoclass:: ProxyPointTarget .. autoclass:: ProxyClusterGeometryData @@ -55,7 +53,6 @@ .. autoclass:: ProxyGenerator .. autoclass:: QBXProxyGenerator -.. autofunction:: partition_by_nodes .. autofunction:: gather_cluster_neighbor_points """ @@ -63,71 +60,6 @@ _DEFAULT_MAX_PARTICLES_IN_BOX = 32 -# {{{ point index partitioning - -def partition_by_nodes( - actx: PyOpenCLArrayContext, places: GeometryCollection, *, - dofdesc: Optional["DOFDescriptorLike"] = None, - tree_kind: Optional[str] = "adaptive-level-restricted", - max_particles_in_box: Optional[int] = None) -> IndexList: - """Generate equally sized ranges of nodes. The partition is created at the - lowest level of granularity, i.e. nodes. This results in balanced ranges - of points, but will split elements across different ranges. - - :arg dofdesc: a :class:`~pytential.symbolic.dof_desc.DOFDescriptor` for - the geometry in *places* which should be partitioned. - :arg tree_kind: if not *None*, it is passed to :class:`boxtree.TreeBuilder`. - :arg max_particles_in_box: value used to control the number of points - in each partition (and thus the number of partitions). See the documentation - in :class:`boxtree.TreeBuilder`. - """ - if dofdesc is None: - dofdesc = places.auto_source - dofdesc = sym.as_dofdesc(dofdesc) - - if max_particles_in_box is None: - max_particles_in_box = _DEFAULT_MAX_PARTICLES_IN_BOX - - lpot_source = places.get_geometry(dofdesc.geometry) - discr = places.get_discretization(dofdesc.geometry, dofdesc.discr_stage) - - if tree_kind is not None: - from pytential.qbx.utils import tree_code_container - tcc = tree_code_container(lpot_source._setup_actx) - - tree, _ = tcc.build_tree()(actx.queue, - particles=flatten( - actx.thaw(discr.nodes()), actx, leaf_class=DOFArray - ), - max_particles_in_box=max_particles_in_box, - kind=tree_kind) - - from boxtree import box_flags_enum - tree = tree.get(actx.queue) - leaf_boxes, = (tree.box_flags & box_flags_enum.HAS_CHILDREN == 0).nonzero() - - indices = np.empty(len(leaf_boxes), dtype=object) - starts = None - - for i, ibox in enumerate(leaf_boxes): - box_start = tree.box_source_starts[ibox] - box_end = box_start + tree.box_source_counts_cumul[ibox] - indices[i] = tree.user_source_ids[box_start:box_end] - else: - if discr.ambient_dim != 2 and discr.dim == 1: - raise ValueError("only curves are supported for 'tree_kind=None'") - - nclusters = max(discr.ndofs // max_particles_in_box, 2) - indices = np.arange(0, discr.ndofs, dtype=np.int64) - starts = np.linspace(0, discr.ndofs, nclusters + 1, dtype=np.int64) - assert starts[-1] == discr.ndofs - - from pytential.linalg import make_index_list - return make_index_list(indices, starts=starts) - -# }}} - - # {{{ proxy points class ProxyPointSource(PointPotentialSource): @@ -178,12 +110,12 @@ class ProxyClusterGeometryData: """ .. attribute:: srcindex - A :class:`~pytential.linalg.IndexList` describing which cluster + A :class:`~pytential.linalg.utils.IndexList` describing which cluster of points each proxy ball was created from. .. attribute:: pxyindex - A :class:`~pytential.linalg.IndexList` describing which proxies + A :class:`~pytential.linalg.utils.IndexList` describing which proxies belong to which cluster. .. attribute:: points diff --git a/pytential/linalg/skeletonization.py b/pytential/linalg/skeletonization.py index ac6fcad4e..ba19afcd5 100644 --- a/pytential/linalg/skeletonization.py +++ b/pytential/linalg/skeletonization.py @@ -22,31 +22,34 @@ from dataclasses import dataclass from typing import ( - Any, Callable, Dict, Hashable, Iterable, Optional, Tuple, Union) + Any, Callable, Dict, Hashable, Iterable, Optional, Tuple, Union, + TYPE_CHECKING) import numpy as np from arraycontext import PyOpenCLArrayContext, Array, ArrayT from pytential import GeometryCollection, sym -from pytential.linalg.utils import IndexList, TargetAndSourceClusterList from pytential.linalg.proxy import ProxyGeneratorBase, ProxyClusterGeometryData +from pytential.linalg.cluster import ClusterTreeLevel, cluster from pytential.linalg.direct_solver_symbolic import ( PROXY_SKELETONIZATION_TARGET, PROXY_SKELETONIZATION_SOURCE, prepare_expr, prepare_proxy_expr) +if TYPE_CHECKING: + from pytential.linalg.utils import IndexList, TargetAndSourceClusterList -__doc__ = """ -.. currentmodule:: pytential.linalg +__doc__ = """ Skeletonization ---------------- +~~~~~~~~~~~~~~~ .. autoclass:: SkeletonizationWrangler .. autoclass:: make_skeletonization_wrangler .. autoclass:: SkeletonizationResult .. autofunction:: skeletonize_by_proxy +.. autofunction:: rec_skeletonize_by_proxy """ @@ -154,8 +157,8 @@ def _apply_weights( actx: PyOpenCLArrayContext, mat: ArrayT, places: GeometryCollection, - tgt_pxy_index: TargetAndSourceClusterList, - cluster_index: IndexList, + tgt_pxy_index: "TargetAndSourceClusterList", + cluster_index: "IndexList", domain: sym.DOFDescriptor) -> ArrayT: """Computes the weights using :func:`_approximate_geometry_waa_magnitude` and multiplies each cluster in *mat* by its corresponding weight. @@ -231,13 +234,13 @@ class SkeletonizationWrangler: A callable that is used to evaluate farfield proxy interactions. This should follow the calling convention of the constructor to - :class:`pytential.symbolic.matrix.P2PClusterMatrixBuilder`. + ``pytential.symbolic.matrix.P2PClusterMatrixBuilder``. .. attribute:: neighbor_cluster_builder A callable that is used to evaluate nearfield neighbour interactions. This should follow the calling convention of the constructor to - :class:`pytential.symbolic.matrix.QBXClusterMatrixBuilder`. + ``pytential.symbolic.matrix.QBXClusterMatrixBuilder``. .. automethod:: evaluate_source_proxy_interaction .. automethod:: evaluate_target_proxy_interaction @@ -292,11 +295,14 @@ def _evaluate_expr(self, # {{{ nearfield - def evaluate_source_neighbor_interaction(self, + def evaluate_source_neighbor_interaction( + self, actx: PyOpenCLArrayContext, places: GeometryCollection, - pxy: ProxyClusterGeometryData, nbrindex: IndexList, *, - ibrow: int, ibcol: int) -> Tuple[np.ndarray, TargetAndSourceClusterList]: + pxy: ProxyClusterGeometryData, nbrindex: "IndexList", *, + ibrow: int, ibcol: int, + ) -> Tuple[np.ndarray, "TargetAndSourceClusterList"]: + from pytential.linalg.utils import TargetAndSourceClusterList nbr_src_index = TargetAndSourceClusterList(nbrindex, pxy.srcindex) eval_mapper_cls = self.neighbor_cluster_builder @@ -307,11 +313,14 @@ def evaluate_source_neighbor_interaction(self, return mat, nbr_src_index - def evaluate_target_neighbor_interaction(self, + def evaluate_target_neighbor_interaction( + self, actx: PyOpenCLArrayContext, places: GeometryCollection, - pxy: ProxyClusterGeometryData, nbrindex: IndexList, *, - ibrow: int, ibcol: int) -> Tuple[np.ndarray, TargetAndSourceClusterList]: + pxy: ProxyClusterGeometryData, nbrindex: "IndexList", *, + ibrow: int, ibcol: int, + ) -> Tuple[np.ndarray, "TargetAndSourceClusterList"]: + from pytential.linalg.utils import TargetAndSourceClusterList tgt_nbr_index = TargetAndSourceClusterList(pxy.srcindex, nbrindex) eval_mapper_cls = self.neighbor_cluster_builder @@ -326,13 +335,17 @@ def evaluate_target_neighbor_interaction(self, # {{{ proxy - def evaluate_source_proxy_interaction(self, + def evaluate_source_proxy_interaction( + self, actx: PyOpenCLArrayContext, places: GeometryCollection, - pxy: ProxyClusterGeometryData, nbrindex: IndexList, *, - ibrow: int, ibcol: int) -> Tuple[np.ndarray, TargetAndSourceClusterList]: - from pytential.collection import add_geometry_to_collection + pxy: ProxyClusterGeometryData, nbrindex: "IndexList", *, + ibrow: int, ibcol: int, + ) -> Tuple[np.ndarray, "TargetAndSourceClusterList"]: + from pytential.linalg.utils import TargetAndSourceClusterList pxy_src_index = TargetAndSourceClusterList(pxy.pxyindex, pxy.srcindex) + + from pytential.collection import add_geometry_to_collection places = add_geometry_to_collection( places, {PROXY_SKELETONIZATION_TARGET: pxy.as_targets()} ) @@ -347,13 +360,17 @@ def evaluate_source_proxy_interaction(self, return mat, pxy_src_index - def evaluate_target_proxy_interaction(self, + def evaluate_target_proxy_interaction( + self, actx: PyOpenCLArrayContext, places: GeometryCollection, - pxy: ProxyClusterGeometryData, nbrindex: IndexList, *, - ibrow: int, ibcol: int) -> Tuple[np.ndarray, TargetAndSourceClusterList]: - from pytential.collection import add_geometry_to_collection + pxy: ProxyClusterGeometryData, nbrindex: "IndexList", *, + ibrow: int, ibcol: int, + ) -> Tuple[np.ndarray, "TargetAndSourceClusterList"]: + from pytential.linalg.utils import TargetAndSourceClusterList tgt_pxy_index = TargetAndSourceClusterList(pxy.srcindex, pxy.pxyindex) + + from pytential.collection import add_geometry_to_collection places = add_geometry_to_collection( places, {PROXY_SKELETONIZATION_SOURCE: pxy.as_sources()} ) @@ -510,10 +527,10 @@ class _ProxyNeighborEvaluationResult: pxy: ProxyClusterGeometryData pxymat: np.ndarray - pxyindex: TargetAndSourceClusterList + pxyindex: "TargetAndSourceClusterList" nbrmat: np.ndarray - nbrindex: TargetAndSourceClusterList + nbrindex: "TargetAndSourceClusterList" def __getitem__(self, i: int) -> Tuple[np.ndarray, np.ndarray]: """ @@ -536,11 +553,11 @@ def _evaluate_proxy_skeletonization_interaction( places: GeometryCollection, proxy_generator: ProxyGeneratorBase, wrangler: SkeletonizationWrangler, - cluster_index: IndexList, *, + cluster_index: "IndexList", *, evaluate_proxy: Callable[..., - Tuple[np.ndarray, TargetAndSourceClusterList]], + Tuple[np.ndarray, "TargetAndSourceClusterList"]], evaluate_neighbor: Callable[..., - Tuple[np.ndarray, TargetAndSourceClusterList]], + Tuple[np.ndarray, "TargetAndSourceClusterList"]], dofdesc: Optional[sym.DOFDescriptor] = None, max_particles_in_box: Optional[int] = None, ) -> _ProxyNeighborEvaluationResult: @@ -551,7 +568,7 @@ def _evaluate_proxy_skeletonization_interaction( if cluster_index.nclusters == 1: raise ValueError("cannot make a proxy skeleton for a single cluster") - from pytential.linalg import gather_cluster_neighbor_points + from pytential.linalg.proxy import gather_cluster_neighbor_points pxy = proxy_generator(actx, dofdesc, cluster_index) nbrindex = gather_cluster_neighbor_points( actx, pxy, @@ -571,7 +588,7 @@ def _skeletonize_block_by_proxy_with_mats( places: GeometryCollection, proxy_generator: ProxyGeneratorBase, wrangler: SkeletonizationWrangler, - tgt_src_index: TargetAndSourceClusterList, *, + tgt_src_index: "TargetAndSourceClusterList", *, id_eps: Optional[float] = None, id_rank: Optional[int] = None, max_particles_in_box: Optional[int] = None ) -> "SkeletonizationResult": @@ -643,7 +660,7 @@ def _skeletonize_block_by_proxy_with_mats( assert R[i].shape == (k, src_mat.shape[1]) assert L[i].shape == (tgt_mat.shape[0], k) - from pytential.linalg import make_index_list + from pytential.linalg import TargetAndSourceClusterList, make_index_list src_skl_index = make_index_list(np.hstack(src_skl_indices), skel_starts) tgt_skl_index = make_index_list(np.hstack(tgt_skl_indices), skel_starts) skel_tgt_src_index = TargetAndSourceClusterList(tgt_skl_index, src_skl_index) @@ -693,13 +710,13 @@ class SkeletonizationResult: .. attribute:: tgt_src_index - A :class:`~pytential.linalg.TargetAndSourceClusterList` representing the - indices in the original matrix :math:`A` that have been skeletonized. + A :class:`~pytential.linalg.utils.TargetAndSourceClusterList` representing + the indices in the original matrix :math:`A` that has been skeletonized. .. attribute:: skel_tgt_src_index - A :class:`~pytential.linalg.TargetAndSourceClusterList` representing a - subset of :attr:`tgt_src_index`, i.e. the skeleton of each cluster of + A :class:`~pytential.linalg.utils.TargetAndSourceClusterList` representing + a subset of :attr:`tgt_src_index`, i.e. the skeleton of each cluster of :math:`A`. These indices can be used to reconstruct the :math:`S` matrix. """ @@ -707,13 +724,13 @@ class SkeletonizationResult: L: np.ndarray R: np.ndarray - tgt_src_index: TargetAndSourceClusterList - skel_tgt_src_index: TargetAndSourceClusterList + tgt_src_index: "TargetAndSourceClusterList" + skel_tgt_src_index: "TargetAndSourceClusterList" # NOTE: these are meant only for testing! They contain the interactions # between the source / target points and their proxies / neighbors. - _src_eval_result: Optional[_ProxyNeighborEvaluationResult] = None - _tgt_eval_result: Optional[_ProxyNeighborEvaluationResult] = None + _src_eval_result: Optional["_ProxyNeighborEvaluationResult"] = None + _tgt_eval_result: Optional["_ProxyNeighborEvaluationResult"] = None def __post_init__(self): if __debug__: @@ -740,7 +757,7 @@ def skeletonize_by_proxy( actx: PyOpenCLArrayContext, places: GeometryCollection, - tgt_src_index: TargetAndSourceClusterList, + tgt_src_index: "TargetAndSourceClusterList", exprs: Union[sym.var, Iterable[sym.var]], input_exprs: Union[sym.var, Iterable[sym.var]], *, domains: Optional[Iterable[Hashable]] = None, @@ -754,7 +771,7 @@ def skeletonize_by_proxy( max_particles_in_box: Optional[int] = None) -> np.ndarray: r"""Evaluate and skeletonize a symbolic expression using proxy-based methods. - :arg tgt_src_index: a :class:`~pytential.linalg.TargetAndSourceClusterList` + :arg tgt_src_index: a :class:`~pytential.linalg.utils.TargetAndSourceClusterList` indicating which indices participate in the skeletonization. :arg exprs: see :func:`make_skeletonization_wrangler`. @@ -762,8 +779,8 @@ def skeletonize_by_proxy( :arg domains: see :func:`make_skeletonization_wrangler`. :arg context: see :func:`make_skeletonization_wrangler`. - :arg approx_nproxy: see :class:`~pytential.linalg.ProxyGenerator`. - :arg proxy_radius_factor: see :class:`~pytential.linalg.ProxyGenerator`. + :arg approx_nproxy: see :class:`~pytential.linalg.proxy.ProxyGenerator`. + :arg proxy_radius_factor: see :class:`~pytential.linalg.proxy.ProxyGenerator`. :arg id_eps: a floating point value used as a tolerance in skeletonizing each block in *tgt_src_index*. @@ -782,14 +799,75 @@ def skeletonize_by_proxy( approx_nproxy=approx_nproxy, radius_factor=proxy_radius_factor) + if wrangler.nrows != 1 or wrangler.ncols != 1: + raise NotImplementedError("support for block matrices") + + from itertools import product skels = np.empty((wrangler.nrows, wrangler.ncols), dtype=object) - for ibrow in range(wrangler.nrows): - for ibcol in range(wrangler.ncols): - skels[ibrow, ibcol] = _skeletonize_block_by_proxy_with_mats( - actx, ibrow, ibcol, places, proxy, wrangler, tgt_src_index, - id_eps=id_eps, id_rank=id_rank, - max_particles_in_box=max_particles_in_box) + for ibrow, ibcol in product(range(wrangler.nrows), range(wrangler.ncols)): + skels[ibrow, ibcol] = _skeletonize_block_by_proxy_with_mats( + actx, ibrow, ibcol, proxy.places, proxy, wrangler, tgt_src_index, + id_eps=id_eps, id_rank=id_rank, + max_particles_in_box=max_particles_in_box) return skels # }}} + + +# {{{ recursive skeletonization by proxy + +def rec_skeletonize_by_proxy( + actx: PyOpenCLArrayContext, + places: GeometryCollection, + + ctree: ClusterTreeLevel, + tgt_src_index: "TargetAndSourceClusterList", + exprs: Union[sym.var, Iterable[sym.var]], + input_exprs: Union[sym.var, Iterable[sym.var]], *, + domains: Optional[Iterable[Hashable]] = None, + context: Optional[Dict[str, Any]] = None, + + approx_nproxy: Optional[int] = None, + proxy_radius_factor: Optional[float] = None, + + id_eps: Optional[float] = None, + id_rank: Optional[int] = None, + max_particles_in_box: Optional[int] = None) -> np.ndarray: + r"""Performs recursive skeletonization based on :func:`skeletonize_by_proxy`. + + :returns: an object :class:`~numpy.ndarray` of :class:`SkeletonizationResult`\ s, + one per level in *ctree*. + """ + + assert ctree.nclusters == tgt_src_index.nclusters + + from pytential.linalg.proxy import QBXProxyGenerator + wrangler = make_skeletonization_wrangler( + places, exprs, input_exprs, + domains=domains, context=context) + proxy = QBXProxyGenerator(places, + approx_nproxy=approx_nproxy, + radius_factor=proxy_radius_factor) + + if wrangler.nrows != 1 or wrangler.ncols != 1: + raise NotImplementedError("support for block matrices") + + from itertools import product + skel_per_level = np.empty(ctree.nlevels, dtype=object) + for i in range(ctree.nlevels): + skels = np.empty((wrangler.nrows, wrangler.ncols), dtype=object) + + for ibrow, ibcol in product(range(wrangler.nrows), range(wrangler.ncols)): + skels[ibrow, ibcol] = _skeletonize_block_by_proxy_with_mats( + actx, ibrow, ibcol, proxy.places, proxy, wrangler, tgt_src_index, + id_eps=id_eps, id_rank=id_rank, + max_particles_in_box=max_particles_in_box) + + skel_per_level[i] = skels + tgt_src_index = cluster(tgt_src_index, ctree) + ctree = ctree.parent() + + return skel_per_level + +# }}} diff --git a/pytential/linalg/utils.py b/pytential/linalg/utils.py index 151b394b4..f6a190c91 100644 --- a/pytential/linalg/utils.py +++ b/pytential/linalg/utils.py @@ -34,16 +34,14 @@ __doc__ = """ -Misc -~~~~ - -.. currentmodule:: pytential.linalg - .. autoclass:: IndexList .. autoclass:: TargetAndSourceClusterList .. autofunction:: make_index_list .. autofunction:: make_index_cluster_cartesian_product +.. autofunction:: make_flat_cluster_diag + +.. autofunction:: interp_decomp """ diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index dfee08d2f..24fa84226 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -245,7 +245,7 @@ def map_call(self, expr): class ClusterMatrixBuilderBase(MatrixBuilderBase): """Evaluate individual clusters of a matrix operator, as defined by a - :class:`~pytential.lingla.TargetAndSourceClusterList`. + :class:`~pytential.linalg.utils.TargetAndSourceClusterList`. Unlike, e.g. :class:`MatrixBuilder`, matrix cluster builders are significantly reduced in scope. They are basically just meant @@ -257,7 +257,8 @@ class ClusterMatrixBuilderBase(MatrixBuilderBase): def __init__(self, actx, dep_expr, other_dep_exprs, dep_source, dep_discr, places, tgt_src_index, context): """ - :arg tgt_src_index: a :class:`~pytential.linalg.TargetAndSourceClusterList` + :arg tgt_src_index: a + :class:`~pytential.linalg.utils.TargetAndSourceClusterList` class describing which clusters are going to be evaluated. """ diff --git a/test/extra_matrix_data.py b/test/extra_matrix_data.py index 62deefd76..eb7ffd4a9 100644 --- a/test/extra_matrix_data.py +++ b/test/extra_matrix_data.py @@ -52,8 +52,8 @@ def get_cluster_index(self, actx, places, dofdesc=None): if max_particles_in_box is None: max_particles_in_box = discr.ndofs // self.approx_cluster_count - from pytential.linalg import partition_by_nodes - cindex = partition_by_nodes(actx, places, + from pytential.linalg.cluster import partition_by_nodes + cindex, ctree = partition_by_nodes(actx, places, dofdesc=dofdesc, tree_kind=self.tree_kind, max_particles_in_box=max_particles_in_box) @@ -74,12 +74,12 @@ def get_cluster_index(self, actx, places, dofdesc=None): from pytential.linalg import make_index_list cindex = make_index_list(subset) - return cindex + return cindex, ctree def get_tgt_src_cluster_index(self, actx, places, dofdesc=None): from pytential.linalg import TargetAndSourceClusterList - cindex = self.get_cluster_index(actx, places, dofdesc=dofdesc) - return TargetAndSourceClusterList(cindex, cindex) + cindex, ctree = self.get_cluster_index(actx, places, dofdesc=dofdesc) + return TargetAndSourceClusterList(cindex, cindex), ctree def get_operator(self, ambient_dim, qbx_forced_limit=_NoArgSentinel): knl = self.knl_class(ambient_dim) diff --git a/test/test_linalg_cluster.py b/test/test_linalg_cluster.py new file mode 100644 index 000000000..b5df0b093 --- /dev/null +++ b/test/test_linalg_cluster.py @@ -0,0 +1,131 @@ +__copyright__ = "Copyright (C) 2022 Alexandru Fikl" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import pytest + +import numpy as np + +from pytential import GeometryCollection +from meshmode.mesh.generation import NArmedStarfish + +from meshmode import _acf # noqa: F401 +from arraycontext import pytest_generate_tests_for_array_contexts +from meshmode.array_context import PytestPyOpenCLArrayContextFactory + +import extra_matrix_data as extra +import logging +logger = logging.getLogger(__name__) + +pytest_generate_tests = pytest_generate_tests_for_array_contexts([ + PytestPyOpenCLArrayContextFactory, + ]) + +CLUSTER_TEST_CASES = [ + extra.CurveTestCase( + name="starfish", + target_order=4, + curve_fn=NArmedStarfish(5, 0.25), + resolutions=[64]), + extra.TorusTestCase( + target_order=4, + resolutions=[1]) + ] + + +# {{{ test_cluster_tree + +@pytest.mark.parametrize(("case", "tree_kind"), [ + (CLUSTER_TEST_CASES[0], None), + (CLUSTER_TEST_CASES[0], "adaptive"), + (CLUSTER_TEST_CASES[0], "adaptive-level-restricted"), + (CLUSTER_TEST_CASES[1], "adaptive"), + ]) +def test_cluster_tree(actx_factory, case, tree_kind, visualize=False): + if visualize: + logging.basicConfig(level=logging.INFO) + + from dataclasses import replace + actx = actx_factory() + case = replace(case, tree_kind=tree_kind) + logger.info("\n%s", case) + + discr = case.get_discretization(actx, case.resolutions[-1], case.target_order) + places = GeometryCollection(discr, auto_where=case.name) + + srcindex, ctree = case.get_cluster_index(actx, places) + assert srcindex.nclusters == ctree.nclusters + + logger.info("nclusters %4d nlevels %4d", srcindex.nclusters, ctree.nlevels) + + if visualize and case.ambient_dim == 2: + import matplotlib.pyplot as plt + fig = plt.figure(figsize=(10, 10), dpi=300) + + from boxtree.visualization import TreePlotter + plotter = TreePlotter(ctree._tree) + plotter.draw_tree(fill=False, edgecolor="black", zorder=10) + plotter.draw_box_numbers() + plotter.set_bounding_box() + + fig.savefig("test_cluster_tree") + + from pytential.linalg.cluster import cluster + for level in range(ctree.nlevels - 1, 0, -1): + logger.info("================ Level %d", ctree.level) + logger.info("partition_box_ids %s", ctree.partition_box_ids) + logger.info("partition_parent_ids %s", ctree.partition_parent_ids) + logger.info("parition_sizes %s", np.diff(srcindex.starts)) + logger.info("partition_parent_map %s", ctree.partition_parent_map) + + next_srcindex = cluster(srcindex, ctree) + next_ctree = ctree.parent() + + assert ctree.level == level + assert srcindex.nclusters == ctree.nclusters + for i, ppm in enumerate(ctree.partition_parent_map): + partition = np.concatenate([srcindex.cluster_indices(j) for j in ppm]) + + assert partition.size == next_srcindex.cluster_size(i) + assert np.allclose(partition, next_srcindex.cluster_indices(i)) + + srcindex = next_srcindex + ctree = next_ctree + + logger.info("================ Level %d", ctree.level) + logger.info("partition_box_ids %s", ctree.partition_box_ids) + logger.info("partition_parent_ids %s", ctree.partition_parent_ids) + logger.info("parition_sizes %s", np.diff(srcindex.starts)) + logger.info("partition_parent_map %s", ctree.partition_parent_map) + assert srcindex.nclusters == ctree.nclusters == 1 + +# }}} + + +if __name__ == "__main__": + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + from pytest import main + main([__file__]) + +# vim: fdm=marker diff --git a/test/test_linalg_proxy.py b/test/test_linalg_proxy.py index 3acb87aea..c67244587 100644 --- a/test/test_linalg_proxy.py +++ b/test/test_linalg_proxy.py @@ -29,7 +29,7 @@ from arraycontext import flatten, unflatten from pytential import bind, sym from pytential import GeometryCollection -from pytential.linalg import ProxyGenerator, QBXProxyGenerator +from pytential.linalg.proxy import ProxyGenerator, QBXProxyGenerator from meshmode.mesh.generation import ellipse, NArmedStarfish from meshmode import _acf # noqa: F401 @@ -211,7 +211,7 @@ def test_partition_points(actx_factory, tree_kind, case, visualize=False): places = GeometryCollection(qbx, auto_where=case.name) density_discr = places.get_discretization(case.name) - mindex = case.get_cluster_index(actx, places) + mindex, _ = case.get_cluster_index(actx, places) expected_indices = np.arange(0, density_discr.ndofs) assert mindex.starts[-1] == density_discr.ndofs @@ -258,7 +258,7 @@ def test_proxy_generator(actx_factory, case, places = GeometryCollection(qbx, auto_where=case.name) density_discr = places.get_discretization(case.name) - cindex = case.get_cluster_index(actx, places) + cindex, _ = case.get_cluster_index(actx, places) generator = proxy_generator_cls(places, approx_nproxy=case.proxy_approx_count, @@ -327,7 +327,7 @@ def test_neighbor_points(actx_factory, case, dofdesc = places.auto_source density_discr = places.get_discretization(dofdesc.geometry) - srcindex = case.get_cluster_index(actx, places) + srcindex, _ = case.get_cluster_index(actx, places) # generate proxy points generator = proxy_generator_cls(places, @@ -336,7 +336,7 @@ def test_neighbor_points(actx_factory, case, pxy = generator(actx, dofdesc, srcindex) # get neighboring points - from pytential.linalg import gather_cluster_neighbor_points + from pytential.linalg.proxy import gather_cluster_neighbor_points nbrindex = gather_cluster_neighbor_points(actx, pxy) pxy = pxy.to_numpy(actx) diff --git a/test/test_linalg_skeletonization.py b/test/test_linalg_skeletonization.py index 33088723d..319fa3773 100644 --- a/test/test_linalg_skeletonization.py +++ b/test/test_linalg_skeletonization.py @@ -115,7 +115,7 @@ def test_skeletonize_symbolic(actx_factory, case, visualize=False): places = GeometryCollection(qbx, auto_where=dd) density_discr = places.get_discretization(dd.geometry, dd.discr_stage) - tgt_src_index = case.get_tgt_src_cluster_index(actx, places, dd) + tgt_src_index, _ = case.get_tgt_src_cluster_index(actx, places, dd) logger.info("nclusters %3d ndofs %7d", tgt_src_index.nclusters, density_discr.ndofs) @@ -124,7 +124,7 @@ def test_skeletonize_symbolic(actx_factory, case, visualize=False): # {{{ wranglers - from pytential.linalg import QBXProxyGenerator + from pytential.linalg.proxy import QBXProxyGenerator from pytential.linalg.skeletonization import make_skeletonization_wrangler proxy_generator = QBXProxyGenerator(places, radius_factor=case.proxy_radius_factor, @@ -157,6 +157,7 @@ def test_skeletonize_symbolic(actx_factory, case, visualize=False): def run_skeletonize_by_proxy(actx, case, resolution, places=None, mat=None, ctol=None, rtol=None, + tgt_src_index=None, suffix="", visualize=False): from pytools import ProcessTimer @@ -166,9 +167,12 @@ def run_skeletonize_by_proxy(actx, case, resolution, if places is None: qbx = case.get_layer_potential(actx, resolution, case.target_order) places = GeometryCollection(qbx, auto_where=dd) + else: + qbx = places.get_geometry(dd.geometry) density_discr = places.get_discretization(dd.geometry, dd.discr_stage) - tgt_src_index = case.get_tgt_src_cluster_index(actx, places, dd) + if tgt_src_index is None: + tgt_src_index, _ = case.get_tgt_src_cluster_index(actx, places, dd) logger.info("nclusters %3d ndofs %7d", tgt_src_index.nclusters, density_discr.ndofs) @@ -185,7 +189,7 @@ def run_skeletonize_by_proxy(actx, case, resolution, logger.info("proxy factor %.2f count %7d", case.proxy_radius_factor, proxy_approx_count) - from pytential.linalg import QBXProxyGenerator + from pytential.linalg.proxy import QBXProxyGenerator from pytential.linalg.skeletonization import make_skeletonization_wrangler proxy_generator = QBXProxyGenerator(places, radius_factor=case.proxy_radius_factor, @@ -234,10 +238,13 @@ def run_skeletonize_by_proxy(actx, case, resolution, logger.info("[time] skeletonization by proxy: %s", p) + def intersect1d(x, y): + return np.where((x.reshape(1, -1) - y.reshape(-1, 1)) == 0)[1] + L, R = skeleton.L, skeleton.R for i in range(tgt_src_index.nclusters): # targets (rows) - bi = np.searchsorted( + bi = intersect1d( tgt_src_index.targets.cluster_indices(i), skeleton.skel_tgt_src_index.targets.cluster_indices(i), ) @@ -247,7 +254,7 @@ def run_skeletonize_by_proxy(actx, case, resolution, tgt_error = la.norm(A - L[i] @ S, ord=2) / la.norm(A, ord=2) # sources (columns) - bj = np.searchsorted( + bj = intersect1d( tgt_src_index.sources.cluster_indices(i), skeleton.skel_tgt_src_index.sources.cluster_indices(i), ) @@ -335,7 +342,7 @@ def run_skeletonize_by_proxy(actx, case, resolution, # }}} - return err_f, (places, mat) + return err_f, (places, mat, skeleton) @pytest.mark.parametrize("case", [ @@ -361,12 +368,28 @@ def test_skeletonize_by_proxy(actx_factory, case, visualize=False): case = replace(case, approx_cluster_count=6, id_eps=1.0e-8) logger.info("\n%s", case) - run_skeletonize_by_proxy( - actx, case, case.resolutions[0], - ctol=6, - # FIXME: why is the 3D error so large? - rtol=10**case.ambient_dim, - visualize=visualize) + dd = sym.DOFDescriptor(case.name, discr_stage=case.skel_discr_stage) + qbx = case.get_layer_potential(actx, case.resolutions[0], case.target_order) + places = GeometryCollection(qbx, auto_where=dd) + + tgt_src_index, ctree = case.get_tgt_src_cluster_index(actx, places, dd) + mat = None + + from pytential.linalg.cluster import cluster + while ctree.level: + logger.info("[%2d/%2d] nclusters %3d", + ctree.level, ctree.nlevels, ctree.nclusters) + + _, (_, mat, skeleton) = run_skeletonize_by_proxy( + actx, case, case.resolutions[0], + ctol=6, + # FIXME: why is the 3D error so large? + rtol=10**case.ambient_dim, + places=places, mat=mat, tgt_src_index=tgt_src_index, + visualize=visualize) + + tgt_src_index = cluster(skeleton.skel_tgt_src_index, ctree) + ctree = ctree.parent() # }}} @@ -430,9 +453,11 @@ def test_skeletonize_by_proxy_convergence( places = mat = None for i in range(id_eps.size): case = replace(case, id_eps=id_eps[i], weighted_proxy=weighted) - rec_error[i], (places, mat) = run_skeletonize_by_proxy( - actx, case, r, places=places, mat=mat, - suffix=f"{suffix}_{i:04d}", visualize=False) + + if not was_zero: + rec_error[i], (places, mat, _) = run_skeletonize_by_proxy( + actx, case, r, places=places, mat=mat, + suffix=f"{suffix}_{i:04d}", visualize=False) was_zero = rec_error[i] == 0.0 eoc.add_data_point(id_eps[i], rec_error[i]) diff --git a/test/test_matrix.py b/test/test_matrix.py index 46fb3dd1e..56e699fcb 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -354,7 +354,7 @@ def test_cluster_builder(actx_factory, ambient_dim, # {{{ matrix - mindex = case.get_tgt_src_cluster_index(actx, places) + mindex, _ = case.get_tgt_src_cluster_index(actx, places) kwargs = dict( dep_expr=sym_u, other_dep_exprs=[], @@ -485,8 +485,8 @@ def test_build_matrix_fixed_stage(actx_factory, logger.info("ndofs: %d", target_discr.ndofs) from pytential.linalg import TargetAndSourceClusterList - itargets = case.get_cluster_index(actx, places, target_dd) - jsources = case.get_cluster_index(actx, places, source_dd) + itargets, _ = case.get_cluster_index(actx, places, target_dd) + jsources, _ = case.get_cluster_index(actx, places, source_dd) mindex = TargetAndSourceClusterList(itargets, jsources) kwargs = dict(