From 32ba92249a5e396bdd8afcaa90a597849851f5e9 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 15 Jul 2020 20:34:12 -0500 Subject: [PATCH 01/14] update proxy tests --- pytential/linalg/proxy.py | 245 ++++++++++++-------------------- test/test_linalg_proxy.py | 285 +++++++++++++++++++------------------- 2 files changed, 235 insertions(+), 295 deletions(-) diff --git a/pytential/linalg/proxy.py b/pytential/linalg/proxy.py index 4183cd0a..534d8aef 100644 --- a/pytential/linalg/proxy.py +++ b/pytential/linalg/proxy.py @@ -28,7 +28,9 @@ import numpy.linalg as la from pytools.obj_array import make_obj_array from pytools import memoize_method, memoize_in + from sumpy.tools import BlockIndexRanges +from boxtree.tools import DeviceDataRecord import loopy as lp from loopy.version import MOST_RECENT_LANGUAGE_VERSION @@ -39,15 +41,26 @@ Proxy Point Generation ~~~~~~~~~~~~~~~~~~~~~~ .. autoclass:: ProxyGenerator -.. autofunction:: partition_by_nodes +.. autofunction:: partition_by_nodes .. autofunction:: gather_block_neighbor_points -.. autofunction:: gather_block_interaction_points + """ # {{{ point index partitioning +def make_block_index(actx, indices, ranges=None): + """Wrap a ``(indices, ranges)`` tuple into a ``BlockIndexRanges``.""" + if ranges is None: + ranges = np.cumsum([0] + [r.size for r in indices]) + indices = np.hstack(indices) + + return BlockIndexRanges(actx.context, + actx.freeze(actx.from_numpy(indices)), + actx.freeze(actx.from_numpy(ranges))) + + def partition_by_nodes(actx, discr, tree_kind="adaptive-level-restricted", max_nodes_in_box=None): """Generate equally sized ranges of nodes. The partition is created at the @@ -82,32 +95,49 @@ def partition_by_nodes(actx, discr, & box_flags_enum.HAS_CHILDREN == 0).nonzero() indices = np.empty(len(leaf_boxes), dtype=np.object) + ranges = 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] - - ranges = actx.from_numpy( - np.cumsum([0] + [box.shape[0] for box in indices]) - ) - indices = actx.from_numpy(np.hstack(indices)) else: - indices = actx.from_numpy(np.arange(0, discr.ndofs, dtype=np.int)) - ranges = actx.from_numpy(np.arange( - 0, - discr.ndofs + 1, - max_nodes_in_box, dtype=np.int)) + indices = np.arange(0, discr.ndofs, dtype=np.int) + ranges = np.arange(0, discr.ndofs + 1, max_nodes_in_box, dtype=np.int) - assert ranges[-1] == discr.ndofs - - return BlockIndexRanges(actx.context, - actx.freeze(indices), actx.freeze(ranges)) + return make_block_index(actx, indices, ranges=ranges) # }}} # {{{ proxy point generator +class BlockProxyPoints(DeviceDataRecord): + """ + .. attribute:: indices + + A :class:`~sumpy.tools.BlockIndexRanges` describing which proxies + belong to which block. + + .. attribute:: points + + A concatenated list of all the proxy points. Can be sliced into + using :attr:`indices` (shape ``(dim, nproxies * nblocks)``). + + .. attribute:: centers + + A list of all the proxy ball centers (shape ``(dim, nblocks)``). + + .. attribute:: radii + + A list of all the proxy ball radii (shape ``(nblocks,)``). + """ + + @property + def nproxy(self): + return self.points[0].shape[0] // self.indices.nblocks + + def _generate_unit_sphere(ambient_dim, approx_npoints): """Generate uniform points on a unit sphere. @@ -212,6 +242,8 @@ class ProxyGenerator(object): def nproxy(self): return self.ref_points.shape[1] + # {{{ proxy ball kernels + @memoize_method def get_kernel(self): if self.radius_factor < 1.0: @@ -281,6 +313,8 @@ class ProxyGenerator(object): return knl + # }}} + def __call__(self, actx, source_dd, indices, **kwargs): """Generate proxy points for each given range of source points in the discretization in *source_dd*. @@ -302,8 +336,7 @@ class ProxyGenerator(object): distance ``pxyradii[i]`` from the range center ``pxycenters[i]``. """ - def _affine_map(v, A, b): - return np.dot(A, v) + b + # {{{ compute proxy centers and radii from pytential import bind, sym source_dd = sym.as_dofdesc(source_dd) @@ -328,33 +361,51 @@ class ProxyGenerator(object): srcranges=indices.ranges, **kwargs) from pytential.utils import flatten_to_numpy - centers = flatten_to_numpy(actx, centers_dev) + centers = np.vstack(flatten_to_numpy(actx, centers_dev)) radii = flatten_to_numpy(actx, radii_dev) - proxies = np.empty(indices.nblocks, dtype=np.object) + + # }}} + + # {{{ build proxy points for each block + + def _affine_map(v, A, b): + return np.dot(A, v) + b + + nproxy = self.nproxy * indices.nblocks + proxies = np.empty((discr.ambient_dim, nproxy), dtype=centers.dtype) + pxy_nr_base = 0 + for i in range(indices.nblocks): - proxies[i] = _affine_map(self.ref_points, + bball = _affine_map(self.ref_points, A=(radii[i] * np.eye(self.ambient_dim)), b=centers[:, i].reshape(-1, 1)) - pxyranges = actx.from_numpy(np.arange( - 0, - proxies.shape[0] * proxies[0].shape[1] + 1, - proxies[0].shape[1], - dtype=indices.ranges.dtype)) + proxies[:, pxy_nr_base:pxy_nr_base + self.nproxy] = bball + pxy_nr_base += self.nproxy + proxies = make_obj_array([ - actx.freeze(actx.from_numpy(np.hstack([p[idim] for p in proxies]))) - for idim in range(self.ambient_dim) + actx.freeze(p) for p in actx.from_numpy(proxies) ]) centers = make_obj_array([ - actx.freeze(centers_dev[idim]) - for idim in range(self.ambient_dim) + actx.freeze(c) for c in centers_dev ]) + pxyindices = np.arange(0, nproxy, dtype=indices.indices.dtype) + pxyranges = np.arange(0, nproxy + 1, self.nproxy) + + assert proxies[0].size == pxyranges[-1] + return BlockProxyPoints( + indices=make_block_index(actx, pxyindices, pxyranges), + points=proxies, + centers=centers, + radii=actx.freeze(radii_dev), + ) + +# }}} - assert pxyranges[-1] == proxies[0].shape[0] - return proxies, actx.freeze(pxyranges), centers, actx.freeze(radii_dev) +# {{{ gather_block_neighbor_points -def gather_block_neighbor_points(actx, discr, indices, pxycenters, pxyradii, +def gather_block_neighbor_points(actx, discr, indices, pxy, max_nodes_in_box=None): """Generate a set of neighboring points for each range of points in *discr*. Neighboring points of a range :math:`i` are defined @@ -363,8 +414,7 @@ def gather_block_neighbor_points(actx, discr, indices, pxycenters, pxyradii, :arg discr: a :class:`meshmode.discretization.Discretization`. :arg indices: a :class:`sumpy.tools.BlockIndexRanges`. - :arg pxycenters: an array containing the center of each proxy ball. - :arg pxyradii: an array containing the radius of each proxy ball. + :arg pxy: a :class:`BlockProxyPoints`. :return: a :class:`sumpy.tools.BlockIndexRanges`. """ @@ -396,17 +446,14 @@ def gather_block_neighbor_points(actx, discr, indices, pxycenters, pxyradii, from boxtree.area_query import AreaQueryBuilder builder = AreaQueryBuilder(actx.context) - query, _ = builder(actx.queue, tree, pxycenters, pxyradii) + query, _ = builder(actx.queue, tree, pxy.centers, pxy.radii) # find nodes inside each proxy ball tree = tree.get(actx.queue) query = query.get(actx.queue) - pxycenters = np.vstack([ - actx.to_numpy(pxycenters[idim]) - for idim in range(discr.ambient_dim) - ]) - pxyradii = actx.to_numpy(pxyradii) + pxycenters = np.vstack([actx.to_numpy(c) for c in pxy.centers]) + pxyradii = actx.to_numpy(pxy.radii) nbrindices = np.empty(indices.nblocks, dtype=np.object) for iproxy in range(indices.nblocks): @@ -420,8 +467,7 @@ def gather_block_neighbor_points(actx, discr, indices, pxycenters, pxyradii, iend = istart + tree.box_source_counts_cumul[iboxes] isources = np.hstack([np.arange(s, e) for s, e in zip(istart, iend)]) - nodes = np.vstack([tree.sources[idim][isources] - for idim in range(discr.ambient_dim)]) + nodes = np.vstack([s[isources] for s in tree.sources]) isources = tree.user_source_ids[isources] # get nodes inside the ball but outside the current range @@ -433,118 +479,7 @@ def gather_block_neighbor_points(actx, discr, indices, pxycenters, pxyradii, nbrindices[iproxy] = indices.indices[isources[mask]] - nbrranges = actx.from_numpy(np.cumsum([0] + [n.shape[0] for n in nbrindices])) - nbrindices = actx.from_numpy(np.hstack(nbrindices)) - - return BlockIndexRanges(actx.context, - actx.freeze(nbrindices), actx.freeze(nbrranges)) - - -def gather_block_interaction_points(actx, places, source_dd, indices, - radius_factor=None, - approx_nproxy=None, - max_nodes_in_box=None): - """Generate sets of interaction points for each given range of indices - in the *source* discretization. For each input range of indices, - the corresponding output range of points is consists of: - - - a set of proxy points (or balls) around the range, which - model farfield interactions. These are constructed using - :class:`ProxyGenerator`. - - - a set of neighboring points that are inside the proxy balls, but - do not belong to the given range, which model nearby interactions. - These are constructed with :func:`gather_block_neighbor_points`. - - :arg places: a :class:`~pytential.symbolic.execution.GeometryCollection`. - :arg source_dd: geometry in *places* for which to generate the - interaction points. This is a - :class:`~pytential.symbolic.primitives.DOFDescriptor` describing - the exact discretization. - :arg indices: a :class:`sumpy.tools.BlockIndexRanges` on the - discretization described by *source_dd*. - - :return: a tuple ``(nodes, ranges)``, where each value is a - :class:`pyopencl.array.Array`. For a range :math:`i`, we can - get the slice using ``nodes[ranges[i]:ranges[i + 1]]``. - """ - - @memoize_in(places, "concat_proxy_and_neighbors") - def knl(): - loopy_knl = lp.make_kernel([ - "{[irange, idim]: 0 <= irange < nranges and \ - 0 <= idim < dim}", - "{[ipxy, ingb]: 0 <= ipxy < npxyblock and \ - 0 <= ingb < nngbblock}" - ], - """ - for irange - <> pxystart = pxyranges[irange] - <> pxyend = pxyranges[irange + 1] - <> npxyblock = pxyend - pxystart - - <> ngbstart = nbrranges[irange] - <> ngbend = nbrranges[irange + 1] - <> nngbblock = ngbend - ngbstart - - <> istart = pxyranges[irange] + nbrranges[irange] - nodes[idim, istart + ipxy] = \ - proxies[idim, pxystart + ipxy] \ - {id_prefix=write_pxy,nosync=write_ngb} - nodes[idim, istart + npxyblock + ingb] = \ - sources[idim, nbrindices[ngbstart + ingb]] \ - {id_prefix=write_ngb,nosync=write_pxy} - ranges[irange + 1] = ranges[irange] + npxyblock + nngbblock - end - """, - [ - lp.GlobalArg("sources", None, - shape=(lpot_source.ambient_dim, "nsources"), dim_tags="sep,C"), - lp.GlobalArg("proxies", None, - shape=(lpot_source.ambient_dim, "nproxies"), dim_tags="sep,C"), - lp.GlobalArg("nbrindices", None, - shape="nnbrindices"), - lp.GlobalArg("nodes", None, - shape=(lpot_source.ambient_dim, "nproxies + nnbrindices")), - lp.ValueArg("nsources", np.int), - lp.ValueArg("nproxies", np.int), - lp.ValueArg("nnbrindices", np.int), - "..." - ], - name="concat_proxy_and_neighbors", - default_offset=lp.auto, - silenced_warnings="write_race(write_*)", - fixed_parameters=dict(dim=lpot_source.ambient_dim), - lang_version=MOST_RECENT_LANGUAGE_VERSION) - - loopy_knl = lp.tag_inames(loopy_knl, "idim*:unr") - loopy_knl = lp.split_iname(loopy_knl, "irange", 128, outer_tag="g.0") - - return loopy_knl - - lpot_source = places.get_geometry(source_dd.geometry) - generator = ProxyGenerator(places, - radius_factor=radius_factor, - approx_nproxy=approx_nproxy) - proxies, pxyranges, pxycenters, pxyradii = \ - generator(actx, source_dd, indices) - - discr = places.get_discretization(source_dd.geometry, source_dd.discr_stage) - neighbors = gather_block_neighbor_points(actx, discr, - indices, pxycenters, pxyradii, - max_nodes_in_box=max_nodes_in_box) - - from meshmode.dof_array import flatten, thaw - ranges = actx.zeros(indices.nblocks + 1, dtype=np.int) - _, (nodes, ranges) = knl()(actx.queue, - sources=flatten(thaw(actx, discr.nodes())), - proxies=proxies, - pxyranges=pxyranges, - nbrindices=neighbors.indices, - nbrranges=neighbors.ranges, - ranges=ranges) - - return actx.freeze(nodes), actx.freeze(ranges) + return make_block_index(actx, nbrindices) # }}} diff --git a/test/test_linalg_proxy.py b/test/test_linalg_proxy.py index 09a1367b..66449f14 100644 --- a/test/test_linalg_proxy.py +++ b/test/test_linalg_proxy.py @@ -151,7 +151,7 @@ def test_proxy_generator(ctx_factory, case, index_sparsity_factor, visualize=Fal case = case.copy(index_sparsity_factor=index_sparsity_factor) logger.info("\n%s", case) - # {{{ generate proxies + # {{{ check proxies qbx = case.get_layer_potential(actx, case.resolutions[-1], case.target_order) places = GeometryCollection(qbx, auto_where=case.name) @@ -160,101 +160,106 @@ def test_proxy_generator(ctx_factory, case, index_sparsity_factor, visualize=Fal srcindices = case.get_block_indices(actx, density_discr, matrix_indices=False) from pytential.linalg.proxy import ProxyGenerator - generator = ProxyGenerator(places) - proxies, pxyranges, pxycenters, pxyradii = \ - generator(actx, places.auto_source, srcindices) + proxies = ProxyGenerator(places)(actx, places.auto_source, srcindices) - proxies = np.vstack([actx.to_numpy(p) for p in proxies]) - pxyranges = actx.to_numpy(pxyranges) - pxycenters = np.vstack([actx.to_numpy(c) for c in pxycenters]) - pxyradii = actx.to_numpy(pxyradii) + from pytential.utils import flatten_to_numpy + proxies = proxies.get(queue) + pxypoints = np.vstack(proxies.points) + pxycenters = np.vstack(proxies.centers) + + srcindices = srcindices.get(queue) + nodes = np.vstack(flatten_to_numpy(actx, density_discr.nodes())) for i in range(srcindices.nblocks): - ipxy = np.s_[pxyranges[i]:pxyranges[i + 1]] + isrc = srcindices.block_indices(i) + ipxy = proxies.indices.block_indices(i) - r = la.norm(proxies[:, ipxy] - pxycenters[:, i].reshape(-1, 1), axis=0) - assert np.allclose(r - pxyradii[i], 0.0, atol=1.0e-14) + r = la.norm(pxypoints[:, ipxy] - pxycenters[:, i].reshape(-1, 1), axis=0) + assert np.allclose(r - proxies.radii[i], 0.0, atol=1.0e-14) + + r = la.norm(nodes[:, isrc] - pxycenters[:, i].reshape(-1, 1), axis=0) + assert np.all(r - proxies.radii[i] < 0.0) # }}} # {{{ visualization - srcindices = srcindices.get(queue) - if visualize: - ambient_dim = places.ambient_dim - if ambient_dim == 2: + if not visualize: + return + + ambient_dim = places.ambient_dim + if ambient_dim == 2: + try: import matplotlib.pyplot as pt + except ImportError: + return + + ci = np.vstack(flatten_to_numpy(actx, + bind(places, sym.expansion_centers(ambient_dim, -1))(actx) + )) + ce = np.vstack(flatten_to_numpy(actx, + bind(places, sym.expansion_centers(ambient_dim, +1))(actx) + )) + r = flatten_to_numpy(actx, + bind(places, sym.expansion_radii(ambient_dim))(actx) + ) + + for i in range(srcindices.nblocks): + isrc = srcindices.block_indices(i) + ipxy = proxies.indices.block_indices(i) + + pt.figure(figsize=(10, 8)) + axis = pt.gca() + for j in isrc: + c = pt.Circle(ci[:, j], r[j], color='k', alpha=0.1) + axis.add_artist(c) + c = pt.Circle(ce[:, j], r[j], color='k', alpha=0.1) + axis.add_artist(c) + + pt.plot(nodes[0], nodes[1], 'ko', ms=2.0, alpha=0.5) + pt.plot(nodes[0, srcindices.indices], nodes[1, srcindices.indices], + 'o', ms=2.0) + pt.plot(nodes[0, isrc], nodes[1, isrc], 'o', ms=2.0) + pt.plot(pxypoints[0, ipxy], pxypoints[1, ipxy], 'o', ms=2.0) + pt.axis("equal") + pt.xlim([-1.5, 1.5]) + pt.ylim([-1.5, 1.5]) + + filename = f"test_proxy_generator_{ambient_dim}d_{i:04}" + pt.savefig(filename, dpi=300) + pt.clf() + else: + from meshmode.discretization.visualization import make_visualizer + from meshmode.mesh.processing import ( # noqa + affine_map, merge_disjoint_meshes) + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + + from meshmode.mesh.generation import generate_icosphere + ref_mesh = generate_icosphere(1, proxies.nproxy) + + # NOTE: this does not plot the actual proxy points + for i in range(srcindices.nblocks): + mesh = affine_map(ref_mesh, + A=(proxies.radii[i] * np.eye(ambient_dim)), + b=pxycenters[:, i].reshape(-1)) + + mesh = merge_disjoint_meshes([mesh, density_discr.mesh]) + discr = Discretization(actx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(10)) - from pytential.utils import flatten_to_numpy - density_nodes = np.vstack(flatten_to_numpy(actx, density_discr.nodes())) - ci = bind(places, sym.expansion_centers(ambient_dim, -1))(actx) - ci = np.vstack(flatten_to_numpy(actx, ci)) - ce = bind(places, sym.expansion_centers(ambient_dim, +1))(actx) - ce = np.vstack(flatten_to_numpy(actx, ce)) - r = bind(places, sym.expansion_radii(ambient_dim))(actx) - r = flatten_to_numpy(actx, r) - - for i in range(srcindices.nblocks): - isrc = srcindices.block_indices(i) - ipxy = np.s_[pxyranges[i]:pxyranges[i + 1]] - - pt.figure(figsize=(10, 8)) - axis = pt.gca() - for j in isrc: - c = pt.Circle(ci[:, j], r[j], color='k', alpha=0.1) - axis.add_artist(c) - c = pt.Circle(ce[:, j], r[j], color='k', alpha=0.1) - axis.add_artist(c) - - pt.plot(density_nodes[0], density_nodes[1], - 'ko', ms=2.0, alpha=0.5) - pt.plot(density_nodes[0, srcindices.indices], - density_nodes[1, srcindices.indices], - 'o', ms=2.0) - pt.plot(density_nodes[0, isrc], density_nodes[1, isrc], - 'o', ms=2.0) - pt.plot(proxies[0, ipxy], proxies[1, ipxy], - 'o', ms=2.0) - pt.xlim([-1.5, 1.5]) - pt.ylim([-1.5, 1.5]) - - filename = "test_proxy_generator_{}d_{:04}.png".format( - ambient_dim, i) - pt.savefig(filename, dpi=300) - pt.clf() - else: - from meshmode.discretization.visualization import make_visualizer - from meshmode.mesh.processing import ( # noqa - affine_map, merge_disjoint_meshes) - from meshmode.discretization import Discretization - from meshmode.discretization.poly_element import \ - InterpolatoryQuadratureSimplexGroupFactory - - from meshmode.mesh.generation import generate_icosphere - ref_mesh = generate_icosphere(1, generator.nproxy) - - # NOTE: this does not plot the actual proxy points - for i in range(srcindices.nblocks): - mesh = affine_map(ref_mesh, - A=(pxyradii[i] * np.eye(ambient_dim)), - b=pxycenters[:, i].reshape(-1)) - - mesh = merge_disjoint_meshes([mesh, density_discr.mesh]) - discr = Discretization(actx, mesh, - InterpolatoryQuadratureSimplexGroupFactory(10)) - - vis = make_visualizer(actx, discr, 10) - filename = "test_proxy_generator_{}d_{:04}.vtu".format( - ambient_dim, i) - vis.write_vtk_file(filename, []) + vis = make_visualizer(actx, discr, 10) + + filename = f"test_proxy_generator_{ambient_dim}d_{i:04}.vtu" + vis.write_vtk_file(filename, []) # }}} @pytest.mark.parametrize("case", PROXY_TEST_CASES) @pytest.mark.parametrize("index_sparsity_factor", [1.0, 0.6]) -def test_interaction_points(ctx_factory, - case, index_sparsity_factor, visualize=False): +def test_neighbor_points(ctx_factory, case, index_sparsity_factor, visualize=False): """Test that neighboring points (inside the proxy balls, but outside the current block/cluster) are actually inside. """ @@ -276,84 +281,84 @@ def test_interaction_points(ctx_factory, # generate proxy points from pytential.linalg.proxy import ProxyGenerator - generator = ProxyGenerator(places) - _, _, pxycenters, pxyradii = generator(actx, places.auto_source, srcindices) + proxies = ProxyGenerator(places)(actx, places.auto_source, srcindices) # get neighboring points - from pytential.linalg.proxy import ( # noqa - gather_block_neighbor_points, - gather_block_interaction_points) - nbrindices = gather_block_neighbor_points(actx, density_discr, - srcindices, pxycenters, pxyradii) - nodes, ranges = gather_block_interaction_points(actx, - places, places.auto_source, srcindices) + from pytential.linalg.proxy import gather_block_neighbor_points + nbrindices = gather_block_neighbor_points(actx, + density_discr, srcindices, proxies) srcindices = srcindices.get(queue) nbrindices = nbrindices.get(queue) + from pytential.utils import flatten_to_numpy + proxies = proxies.get(queue) + pxypoints = np.vstack(proxies.points) + pxycenters = np.vstack(proxies.centers) + nodes = np.vstack(flatten_to_numpy(actx, density_discr.nodes())) + for i in range(srcindices.nblocks): isrc = srcindices.block_indices(i) inbr = nbrindices.block_indices(i) assert not np.any(np.isin(inbr, isrc)) + r = la.norm(nodes[:, inbr] - pxycenters[:, i].reshape(-1, 1), axis=0) + assert np.all(r - proxies.radii[i] < 0.0) + # }}} # {{{ visualize - from pytential.utils import flatten_to_numpy - if visualize: - ambient_dim = places.ambient_dim - if ambient_dim == 2: + if not visualize: + return + + ambient_dim = places.ambient_dim + if ambient_dim == 2: + try: import matplotlib.pyplot as pt - density_nodes = flatten_to_numpy(actx, density_discr.nodes()) - nodes = flatten_to_numpy(actx, nodes) - ranges = actx.to_numpy(ranges) - - for i in range(srcindices.nblocks): - isrc = srcindices.block_indices(i) - inbr = nbrindices.block_indices(i) - iall = np.s_[ranges[i]:ranges[i + 1]] - - pt.figure(figsize=(10, 8)) - pt.plot(density_nodes[0], density_nodes[1], - 'ko', ms=2.0, alpha=0.5) - pt.plot(density_nodes[0][srcindices.indices], - density_nodes[1][srcindices.indices], - 'o', ms=2.0) - pt.plot(density_nodes[0][isrc], density_nodes[1][isrc], - 'o', ms=2.0) - pt.plot(density_nodes[0][inbr], density_nodes[1][inbr], - 'o', ms=2.0) - pt.plot(nodes[0][iall], nodes[1][iall], - 'x', ms=2.0) - pt.xlim([-1.5, 1.5]) - pt.ylim([-1.5, 1.5]) - - filename = "test_area_query_{}d_{:04}.png".format(ambient_dim, i) - pt.savefig(filename, dpi=300) - pt.clf() - elif ambient_dim == 3: - from meshmode.discretization.visualization import make_visualizer - marker = np.empty(density_discr.ndofs) - - for i in range(srcindices.nblocks): - isrc = srcindices.block_indices(i) - inbr = nbrindices.block_indices(i) - - marker.fill(0.0) - marker[srcindices.indices] = 0.0 - marker[isrc] = -42.0 - marker[inbr] = +42.0 - - from meshmode.dof_array import unflatten - marker_dev = unflatten(actx, density_discr, actx.from_numpy(marker)) - - vis = make_visualizer(actx, density_discr, 10) - filename = "test_area_query_{}d_{:04}.vtu".format(ambient_dim, i) - vis.write_vtk_file(filename, [ - ("marker", marker_dev), - ]) + except ImportError: + return + + for i in range(srcindices.nblocks): + isrc = srcindices.block_indices(i) + inbr = nbrindices.block_indices(i) + + pt.figure(figsize=(10, 10)) + pt.plot(nodes[0], nodes[1], 'ko', ms=2.0, alpha=0.5) + pt.plot(nodes[0, srcindices.indices], nodes[1, srcindices.indices], + 'o', ms=2.0) + pt.plot(nodes[0, isrc], nodes[1, isrc], 'o', ms=2.0) + pt.plot(nodes[0, inbr], nodes[1, inbr], 'o', ms=2.0) + pt.axis("equal") + pt.xlim([-1.5, 1.5]) + pt.ylim([-1.5, 1.5]) + + filename = f"test_area_query_{ambient_dim}d_{i:04}" + pt.savefig(filename, dpi=300) + pt.clf() + elif ambient_dim == 3: + from meshmode.discretization.visualization import make_visualizer + marker = np.empty(density_discr.ndofs) + + for i in range(srcindices.nblocks): + isrc = srcindices.block_indices(i) + inbr = nbrindices.block_indices(i) + + marker.fill(0.0) + marker[srcindices.indices] = 0.0 + marker[isrc] = -42.0 + marker[inbr] = +42.0 + + from meshmode.dof_array import unflatten + marker_dev = unflatten(actx, density_discr, actx.from_numpy(marker)) + + vis = make_visualizer(actx, density_discr, 10) + filename = "test_area_query_{}d_{:04}.vtu".format(ambient_dim, i) + vis.write_vtk_file(filename, [ + ("marker", marker_dev), + ]) + # }}} -- GitLab From d5cc26b05b21c9da18a32f20c6e4f5d7b9c86701 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 15 Jul 2020 20:34:35 -0500 Subject: [PATCH 02/14] add first stab at skeletonization --- pytential/linalg/skeletonization.py | 413 ++++++++++++++++++++++++++++ 1 file changed, 413 insertions(+) create mode 100644 pytential/linalg/skeletonization.py diff --git a/pytential/linalg/skeletonization.py b/pytential/linalg/skeletonization.py new file mode 100644 index 00000000..06847ed4 --- /dev/null +++ b/pytential/linalg/skeletonization.py @@ -0,0 +1,413 @@ +from __future__ import division, absolute_import + +__copyright__ = "Copyright (C) 2018-2020 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 contextlib import contextmanager + +import numpy as np +import numpy.linalg as la + +import pyopencl as cl + +from pytools.obj_array import make_obj_array +from pytools import Record + +from sumpy.tools import MatrixBlockIndexRanges, BlockIndexRanges + + +# {{{ helpers + +def interp_decomp(A, rank, eps): + """Wrapper for :func:`~scipy.linalg.interpolative.interp_decomp` that + always has the same output signature. + + :return: a tuple ``(k, idx, interp)`` containing the numerical rank, + the column indices and the resulting interpolation matrix. + """ + + import scipy.linalg.interpolative as sli # pylint:disable=no-name-in-module + if rank is None: + k, idx, proj = sli.interp_decomp(A, eps) + else: + idx, proj = sli.interp_decomp(A, rank) + k = rank + + # NOTE: fix should be in scipy 1.2.0 + # https://github.com/scipy/scipy/pull/9125 + if k == A.shape[1]: + proj = np.empty((k, 0), dtype=proj.dtype) + + interp = sli.reconstruct_interp_matrix(idx, proj) + return k, idx, interp + + +def make_block_diag(blk, blkindices): + """Construct a block diagonal matrix from a linear representation of + the matrix blocks. + """ + + nblocks = blkindices.nblocks + diag = np.full((nblocks, nblocks), 0, dtype=np.object) + + for i in range(nblocks): + diag[i, i] = blkindices.block_take(blk, i) + + return diag + + +class SkeletonizedBlock(Record): + """ + .. attribute:: L + .. attribute:: R + .. attribute:: sklindices + """ + +# }}} + + +# {{{ evaluation + +class QBXForcedLimitReplacer(IdentityMapper): + def __init__(self, qbx_forced_limit=None): + super(QBXForcedLimitReplacer, self).__init__() + self.qbx_forced_limit = qbx_forced_limit + + def map_int_g(self, expr): + return expr.copy(qbx_forced_limit=self.qbx_forced_limit) + + +class DOFDescReplacer(ToTargetTagger): + def _default_dofdesc(self, dofdesc): + return self.default_where + + +def prepare_expr(places, expr, auto_where=None, qbx_forced_limit=None): + from pytential.symbolic.execution import _prepare_auto_where + auto_where = _prepare_auto_where(auto_where, places=places) + + from pytential.symbolic.execution import _prepare_expr + expr = _prepare_expr(places, expr, auto_where=auto_where) + + expr = QBXForcedLimitReplacer(qbx_forced_limit)(expr) + expr = DOFDescReplacer(auto_where[0], auto_where[1])(expr) + + return expr + + +class BlockEvaluationWrangler: + """ + .. automethod:: evaluate_source_farfield + .. automethod:: evaluate_target_farfield + .. automethod:: evaluate_nearfield + """ + + def __init__(self, exprs, input_exprs, domains, + context=None, + weighted_farfield=None, + farfield_block_builder=None, + nearfield_block_builder=None): + """ + :arg exprs: an :class:`~numpy.ndarray` of expressions (layer potentials) + that correspond to the output blocks of the matrix. + :arg input_exprs: a :class:`list` of densities that correspond to the + input blocks of the matrix. + :arg domains: a :class:`list` of the same length as *input_exprs* + defining the domain of each input. + + :arg context: a :class:`dict` with additional parameters required to + evaluated the expressions. + + :arg weighted_farfield: a :class:`tuple` containing two :class:`bool`\ s + which turn on the weighing proxy interactions by the corresponding + quadrature weights. + :arg farfield_block_builder: + a :class:`pytential.symbolic.matrix.MatrixBlockBuilderBase` that + is used to evaluate farfield proxy interactions. + :arg nearfield_block_builder: + a :class:`pytential.symbolic.matrix.MatrixBlockBuilderBase` that + is used to evaluate nearfield neighbour interactions. + """ + if context is None: + context = {} + + self.exprs = exprs + self.input_exprs = input_exprs + self.domains = domains + self.context = context + + if weighted_farfield is None: + weighted_source = True + weighted_target = False + elif isinstance(weighted_farfield, bool): + weighted_source = weighted_target = weighted_farfield + elif isinstance(weighted_farfield, (list, tuple)): + weighted_source, weighted_target = weighted_farfield + else: + raise ValueError("unknown value for weighting: `{}`".format( + weighted_farfield)) + + self.weighted_source = weighted_source + self.weighted_target = weighted_target + + self.nearfield_block_builder = nearfield_block_builder + if self.nearfield_block_builder is None: + from pytential.symbolic.matrix import NearFieldBlockBuilder + self.nearfield_block_builder = NearFieldBlockBuilder + + self.farfield_block_builder = farfield_block_builder + if self.farfield_block_builder is None: + from pytential.symbolic.matrix import FarFieldBlockBuilder + self.farfield_block_builder = FarFieldBlockBuilder + + def _evaluate(self, actx, places, builder_cls, + expr, idomain, index_set, auto_where, **kwargs): + domain = self.domains[idomain] + dep_source = places.get_geometry(domain.geometry) + dep_discr = places.get_discretization(domain.geometry, domain.discr_stage) + + expr = prepare_expr(places, expr, auto_where=auto_where) + builder = builder_cls(actx, + dep_expr=self.input_exprs[idomain], + other_dep_exprs=( + self.input_exprs[:idomain] + + self.input_exprs[idomain+1:]), + dep_source=dep_source, + dep_discr=dep_discr, + places=places, + index_set=index_set, + context=self.context, + **kwargs) + + return builder(expr) + + def evaluate_source_farfield(self, + actx, places, ibrow, ibcol, index_set, auto_where=None): + return self._evaluate(actx, places, + self.farfield_block_builder, + expr, ibcol, index_set, auto_where, + weighted=self.weighted_source, + exclude_self=False) + + def evaluate_target_farfield(self, + actx, places, ibrow, ibcol, index_set, auto_where=None): + return self._evaluate(actx, places, + self.farfield_block_builder, + expr, ibcol, index_set, auto_where, + weighted=self.weighted_target, + exclude_self=False) + + def evaluate_nearfield(self, + actx, places, ibrow, ibcol, index_set, auto_where=None): + return self._evaluate(actx, places, + self.nearfield_block_builder, + self.exprs[ibrow], ibcol, index_set, auto_where) + +# }}} + + +# {{{ skeletonize_by_proxy + +@contextmanager +def add_to_geometry_collection(places, proxy): + # NOTE: this is a bit of a hack to keep all the caches in `places` and + # just add the proxy points to it, since otherwise all the DISCRETIZATION + # scope stuff would be recomputed all the time + + try: + previous_cse_cache = places._get_cache("cse") + places.places["proxy"] = proxy + yield places + finally: + del places.places["proxy"] + # NOTE: this is meant to make sure that proxy-related things don't + # get cached over multiple runs and lead to some explosion of some sort + places.cache["cse"] = previous_cse_cache + + +def make_block_proxy_skeleton(actx, places, proxy, wrangler, indices, + ibrow, ibcol, source_or_target, + max_particles_in_box=None): + """Builds a block matrix that can be used to skeletonize the rows + (targets) of the symbolic matrix block described by ``(ibrow, ibcol)``. + """ + + if source_or_target == "source": + from pytential.target import PointsTarget as ProxyPoints + + swap_arg_order = lambda x, y: (x, y) + evaluate_farfield = wrangler.evaluate_source_farfield + elif source_or_target == "target": + from pytential.source import PointPotentialSource as ProxyPoints + + swap_arg_order = lambda x, y: (y, x) + evaluate_farfield = wrangler.evaluate_target_farfield + else: + raise ValueError(f"unknown value: '{source_or_target}'") + + # {{{ generate proxies + + domain = wrangler.domains[ibcol] + dep_source = places.get_geometry(domain.geometry) + dep_discr = places.get_discretization(domain.geometry, domain.discr_stage) + + pxy = proxy(actx, domain, indices) + + # }}} + + # {{{ evaluate (farfield) proxy interactions + + with add_to_geometry_collection(places, ProxyPoints(pxy.points)): + pxyindices = MatrixBlockIndexRanges(actx.context, + *swap_arg_order(pxy.indices, indices)) + pxymat = evaluate_farfield(actx, pxyplaces, ibrow, ibcol, pxyindices, + auto_where=swap_arg_order(domain, "proxy")) + + if indices.nblocks == 1: + # TODO: evaluate nearfield at the root level? + return make_block_diag(pxymat, pxyindices.get(actx.queue)) + + # }}} + + # {{{ evaluate (nearfield) neighbor interactions + + nbrindices = gather_block_neighbor_points( + actx, dep_discr, indices, pxy, + max_particles_in_box=max_particles_in_box) + + nbrindices = MatrixBlockIndexRanges(actx.context, + *swap_arg_order(nbrindices, indices)) + nbrmat = wrangler.evaluate_nearfield(actx, places, ibrow, ibcol, nbrindices) + + # }}} + + # {{{ concatenate everything to get the blocks ready for ID-ing + + pxyindices = pxyindices.get(queue) + nbrindices = nbrindices.get(queue) + + pxyblk = np.full((indices.nblocks, indices.nblocks), 0, dtype=np.object) + for i in range(indices.nblocks): + pxyblk[i, i] = np.hstack(swap_arg_order( + pxyindices.block_take(pxymat, i) + nbrindices.block_take(nbrmat, i), + )) + # }}} + + return pxyblk + + +def skeletonize_block_by_proxy(actx, + ibcol, ibrow, places, proxy, wrangler, blkindices, + id_eps=None, id_rank=None, + tree_max_particles_in_box=None): + r""" + :arg places: a :class:`~meshmode.array_context.ArrayContext`. + :arg proxy: a :class:`~pytential.linalg.proxy.ProxyGenerator`. + :arg wrangler: a :class:`BlockEvaluationWrangler`. + :arg blkindices: a :class:`~sumpy.tools.MatrixBlockIndexRanges`. + + :returns: a tuple ``(L, R, sklindices)`` encoding the block-by-block + decompression of the matrix represented *wrangler*. :math:`L` and + :math:`R` are :math:`n \times n` diagonal block matrix, where + :math:`n` is :attr:`~sumpy.tools.MatrixBlockIndexRanges.nblocks``. The + ``sklindices`` array contains the remaining (skeleton) nodes from + ``blkindices`` after compression. + """ + + L = np.full((blkindices.nblocks, blkindices.nblocks), 0, dtype=np.object) + R = np.full((blkindices.nblocks, blkindices.nblocks), 0, dtype=np.object) + + if blkindices.nblocks == 1: + L[0, 0] = np.eye(blkindices.row.indices.size) + R[0, 0] = np.eye(blkindices.col.indices.size) + + return L, R, blkindices + + # construct proxy matrices to skeletonize + src_mat = build_source_skeleton_matrix(queue, + places, proxy, wrangler, blkindices.col, 0, 0, + max_particles_in_box=tree_max_particles_in_box) + tgt_mat = build_target_skeleton_matrix(queue, + places, proxy, wrangler, blkindices.row, 0, 0, + max_particles_in_box=tree_max_particles_in_box) + + src_skl_indices = np.empty(blkindices.nblocks, dtype=np.object) + tgt_skl_indices = np.empty(blkindices.nblocks, dtype=np.object) + skl_ranges = np.zeros(blkindices.nblocks + 1, dtype=np.int) + + src_indices = blkindices.col.get(queue) + tgt_indices = blkindices.row.get(queue) + + for i in range(blkindices.nblocks): + k = id_rank + + assert not np.any(np.isnan(src_mat[i, i])), "block {}".format(i) + assert not np.any(np.isinf(src_mat[i, i])), "block {}".format(i) + assert not np.any(np.isnan(tgt_mat[i, i])), "block {}".format(i) + assert not np.any(np.isinf(tgt_mat[i, i])), "block {}".format(i) + + # skeletonize target points + k, idx, interp = _interp_decomp(tgt_mat[i, i].T, k, id_eps) + assert k > 0 + + L[i, i] = interp.T + tgt_skl_indices[i] = tgt_indices.block_indices(i)[idx[:k]] + + # skeletonize source points + k, idx, interp = _interp_decomp(src_mat[i, i], k, id_eps) + assert k > 0 + + R[i, i] = interp + src_skl_indices[i] = src_indices.block_indices(i)[idx[:k]] + + skl_ranges[i + 1] = skl_ranges[i] + k + assert R[i, i].shape == (k, src_mat[i, i].shape[1]) + assert L[i, i].shape == (tgt_mat[i, i].shape[0], k) + + src_skl_indices = _to_block_index(queue, src_skl_indices, skl_ranges) + tgt_skl_indices = _to_block_index(queue, tgt_skl_indices, skl_ranges) + skl_indices = MatrixBlockIndexRanges(queue.context, + tgt_skl_indices, + src_skl_indices) + + return L, R, skl_indices + + +def skeletonize_by_proxy(actx, places, proxy, wrangler, blkindices, + id_eps=None, id_rank=None, + max_particles_in_box=None): + + nrows = len(wrangler.exprs) + ncols = len(wrangler.input_exprs) + skel = np.empty((nrows, ncols), dtype=np.object) + + for ibrow in range(nrows): + for ibcol in range(ncols) + skel[ibrow, ibcol] = skeletonize_block_by_proxy(actx, + ibrow, ibcol, proxy, wrangler, blkindices, + id_eps=id_eps, id_rank=id_rank, + max_particles_in_box=max_particles_in_box) + + return skel + +# }}} -- GitLab From 7a1af1c2b2c7169d9a5f9bec7a9da466d4259959 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 15 Jul 2020 21:25:35 -0500 Subject: [PATCH 03/14] add a test case for skeletonization --- pytential/linalg/skeletonization.py | 35 ++++++++- test/extra_matrix_data.py | 91 +++++++++++++++++++++- test/test_linalg_proxy.py | 117 +++++++++++++++++++++++++++- test/test_matrix.py | 3 +- 4 files changed, 239 insertions(+), 7 deletions(-) diff --git a/pytential/linalg/skeletonization.py b/pytential/linalg/skeletonization.py index 06847ed4..95efa15d 100644 --- a/pytential/linalg/skeletonization.py +++ b/pytential/linalg/skeletonization.py @@ -185,7 +185,6 @@ class BlockEvaluationWrangler: dep_source = places.get_geometry(domain.geometry) dep_discr = places.get_discretization(domain.geometry, domain.discr_stage) - expr = prepare_expr(places, expr, auto_where=auto_where) builder = builder_cls(actx, dep_expr=self.input_exprs[idomain], other_dep_exprs=( @@ -202,6 +201,7 @@ class BlockEvaluationWrangler: def evaluate_source_farfield(self, actx, places, ibrow, ibcol, index_set, auto_where=None): + expr = prepare_expr(places, expr, auto_where=auto_where) return self._evaluate(actx, places, self.farfield_block_builder, expr, ibcol, index_set, auto_where, @@ -210,6 +210,7 @@ class BlockEvaluationWrangler: def evaluate_target_farfield(self, actx, places, ibrow, ibcol, index_set, auto_where=None): + expr = prepare_expr(places, expr, auto_where=auto_where) return self._evaluate(actx, places, self.farfield_block_builder, expr, ibcol, index_set, auto_where, @@ -222,6 +223,38 @@ class BlockEvaluationWrangler: self.nearfield_block_builder, self.exprs[ibrow], ibcol, index_set, auto_where) + +def make_block_evaluation_wrangler(places, exprs, input_exprs, + domains=None, context=None, + _weighted_farfield=None, + _farfield_block_builder=None, + _nearfield_block_builder=None): + + if not is_obj_array(exprs): + exprs = make_obj_array([exprs]) + + try: + input_exprs = list(input_exprs) + except TypeError: + input_exprs = [input_exprs] + + from pytential.symbolic.execution import _prepare_auto_where + auto_where = _prepare_auto_where(auto_where, places) + from pytential.symbolic.execution import _prepare_domains + domains = _prepare_domains(len(input_exprs), places, domains, auto_where[0]) + + if context is None: + context = {} + + return BlockEvaluationWrangler( + exprs=exprs, + input_exprs=input_exprs, + domains=domains, + context=context, + weighted_farfield=_weighted_farfield, + farfield_block_builder=_farfield_block_builder, + nearfield_block_builder=_nearfield_block_builder) + # }}} diff --git a/test/extra_matrix_data.py b/test/extra_matrix_data.py index 983c815e..bd78f8a9 100644 --- a/test/extra_matrix_data.py +++ b/test/extra_matrix_data.py @@ -21,23 +21,106 @@ def max_block_error(mat, blk, index_set, p=None): return error + +def skeletonization_error(mat, skel, srcindices, p=None): + r"""Matrix errors are computed as follows: + + .. math:: + + \begin{aligned} + \epsilon_{l, ij} = \max_{i \ne j} \|A_{ij} - L_{ii} S_{ij}\|_p, \\ + \epsilon_{r, ij} = \max_{i \ne j} \|A_{ij} - S_{ij} R_{jj}\|_p, \\ + \epsilon_f = \|A - L S R\|_p. + \end{aligned} + + :arg mat: dense matrix. + :arg skel: a :class:`~pytential.linalg.skeletonization.SkeletonizedBlock`. + :arg p: norm type (follows :func:`~numpy.linalg.norm` rules for matrices). + + :returns: a tuple ``(err_l, err_r, err_f)`` of the left, right and full + matrix errors. + """ + from itertools import product + + L = skel.L + R = skel.R + sklindices = skel.sklindices + + def mnorm(x, y): + return la.norm(x - y, ord=p) / la.norm(x, ord=p) + + # build block matrices + nblocks = srcindices.nblocks + S = np.full((nblocks, nblocks), 0, dtype=np.object) + A = np.full((nblocks, nblocks), 0, dtype=np.object) + + def block_indices(blk, i): + return blk.indices[blk.ranges[i]:blk.ranges[i + 1]] + + # compute max block error + err_l = np.empty((nblocks, nblocks)) + err_l = np.empty((nblocks, nblocks)) + for i, j in product(range(nblocks), repeat=2): + if i == j: + continue + + # full matrix indices + f_tgt = block_indices(srcindices.row, i) + f_src = block_indices(srcindices.col, j) + # skeleton matrix indices + s_tgt = block_indices(sklindices.row, i) + s_src = block_indices(sklindices.col, j) + + S[i, j] = mat[np.ix_(s_tgt, s_src)] + A[i, j] = mat[np.ix_(f_tgt, f_src)] + + blk = mat[np.ix_(s_tgt, f_src)] + err_l[i, j] = mnorm(A[i, j], L[i, i].dot(blk)) + + blk = mat[np.ix_(f_tgt, s_src)] + err_r[i, j] = mnorm(A[i, j], blk.dot(R[j, j])) + + # compute full matrix error + from pytential.symbolic.execution import _bmat + A = _bmat(A, dtype=mat.dtype) + L = _bmat(L, dtype=mat.dtype) + S = _bmat(S, dtype=mat.dtype) + R = _bmat(R, dtype=mat.dtype) + + assert L.shape == (A.shape[0], S.shape[0]) + assert R.shape == (S.shape[1], A.shape[1]) + err_f = mnorm(A, L.dot(S.dot(R))) + + return err_l, err_r, err_f + # }}} # {{{ MatrixTestCase class MatrixTestCaseMixin: + # operators + op_type = "scalar" + # disable fmm for matrix tests + fmm_backend = None + # partitioning approx_block_count = 10 max_particles_in_box = None tree_kind = "adaptive-level-restricted" index_sparsity_factor = 1.0 - # operators - op_type = "scalar" + # proxy + proxy_radius_factor = None + proxy_approx_count = None - # disable fmm for matrix tests - fmm_backend = None + # skeletonization + id_eps = 1.0e-8 + skel_discr_stage = sym.QBX_SOURCE_STAGE2 + + weighted_farfield = None + farfield_block_builder = None + nearfield_block_builder = None def get_block_indices(self, actx, discr, matrix_indices=True): max_particles_in_box = self.max_particles_in_box diff --git a/test/test_linalg_proxy.py b/test/test_linalg_proxy.py index 66449f14..fa460141 100644 --- a/test/test_linalg_proxy.py +++ b/test/test_linalg_proxy.py @@ -46,6 +46,8 @@ logger = logging.getLogger(__name__) logging.basicConfig(level=logging.INFO) +# {{{ plot_partition_indices + def plot_partition_indices(actx, discr, indices, **kwargs): try: import matplotlib.pyplot as pt @@ -54,7 +56,7 @@ def plot_partition_indices(actx, discr, indices, **kwargs): indices = indices.get(actx.queue) args = [ - kwargs.get("tree_kind", "linear").replace("-", "_"), + str(kwargs.get("tree_kind", "linear")).replace("-", "_"), kwargs.get("discr_stage", "stage1"), discr.ambient_dim ] @@ -98,6 +100,8 @@ def plot_partition_indices(actx, discr, indices, **kwargs): ("marker", marker) ]) +# }}} + PROXY_TEST_CASES = [ extra.CurveTestCase( @@ -110,6 +114,8 @@ PROXY_TEST_CASES = [ ] +# {{{ test_partition_points + @pytest.mark.skip(reason="only useful with visualize=True") @pytest.mark.parametrize("tree_kind", ['adaptive', None]) @pytest.mark.parametrize("case", PROXY_TEST_CASES) @@ -136,6 +142,10 @@ def test_partition_points(ctx_factory, tree_kind, case, visualize=False): # }}} +# }}} + + +# {{{ test_proxy_generator @pytest.mark.parametrize("case", PROXY_TEST_CASES) @pytest.mark.parametrize("index_sparsity_factor", [1.0, 0.6]) @@ -256,6 +266,10 @@ def test_proxy_generator(ctx_factory, case, index_sparsity_factor, visualize=Fal # }}} +# }}} + + +# {{{ test_neighbor_points @pytest.mark.parametrize("case", PROXY_TEST_CASES) @pytest.mark.parametrize("index_sparsity_factor", [1.0, 0.6]) @@ -361,6 +375,107 @@ def test_neighbor_points(ctx_factory, case, index_sparsity_factor, visualize=Fal # }}} +# }}} + + +# {{{ test_skeletonize_by_proxy + +@pytest.mark.parametrize("case", PROXY_TEST_CASES) +def test_skeletonize_by_proxy(ctx_factory, case, visualize=False): + """Test single-level level skeletonization accuracy.""" + + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + actx = PyOpenCLArrayContext(queue) + + case = case.copy(nblocks=6, op_type="scalar") + logger.info("\n%s", case) + + # {{{ geometry + + dd = sym.DOFDescriptor(case.name, discr_stage=case.skel_discr_stage) + qbx = case.get_layer_potential(ctx, case.resolutions[0], case.target_order) + places = GeometryCollection(qbx, auto_where=dd) + + density_discr = places.get_discretization(dd.geometry, dd.discr_stage) + srcindices = case.get_block_indices(actx, density_discr) + + logger.info("nblocks %3d ndofs %7d", srcindices.nblocks, density_discr.ndofs) + + # }}} + + # {{{ wranglers + + from pytential.linalg.proxy import ProxyGenerator + from pytential.linalg.skeletonization import make_block_evaluation_wrangler + proxy = ProxyGenerator(places, + radius_factor=case.proxy_radius_factor, + approx_nproxy=case.proxy_approx_count) + + sym_u, sym_op = case.get_operator(places.ambient_dim) + wrangler = make_block_evaluation_wrangler(places, sym_op, sym_u, + domains=None, + context=case.knl_concrete_kwargs, + _weighted_farfield=case.weighted_farfield, + _farfield_block_builder=case.farfield_block_builder, + _nearfield_block_builder=case.nearfield_block_builder) + + # }}} + + # {{{ check skeletonize + + # dense matrix + from pytential.symbolic.execution import build_matrix + mat = actx.to_numpy( + build_matrix(actx, places, sym_op, sym_u, + context=case.knl_concrete_kwargs) + ) + + # skeleton + from pytential.linalg.skeletonize import skeletonize_by_proxy + skeleton = skeletonize_by_proxy(actx, + places, proxy, wrangler, srcindices, + id_eps=case.id_eps, + max_particles_in_box=case.max_particles_in_box) + + err_l, err_r, err_f = extra.skeletonization_error(mat, skeleton, srcindices) + + # }}} + + # {{{ visualize + + if not visualize + return + + def plot_skeleton(isrc, iskl, name): + pt.figure(figsize=(10, 10), dpi=300) + pt.plot(sources[0][isrc.indices], sources[1][isrc.indices], + 'ko', alpha=0.5) + + for i in range(blkindices.nblocks): + iblk = iskl.block_indices(i) + pt.plot(sources[0][iblk], sources[1][iblk], 'o') + + pt.savefig(f"test_skeletonize_by_proxy_{name}") + pt.clf() + + if places.ambient_dim == 2: + sources = flatten_to_numpy(actx, density_discr.nodes()) + srcindices = srcindices.get(queue) + sklindices = sklindices.get(queue) + + plot_skeleton(srcindices.col, sklindices.col, "sources") + plot_skeleton(srcindices.row, sklindices.row, "targets") + else: + # TODO: would be nice to figure out a way to visualize some of these + # skeletonization results in 3D. Probably need to teach the + # visualizers to spit out point clouds + pass + + # }}} + +# }}} + if __name__ == "__main__": import sys diff --git a/test/test_matrix.py b/test/test_matrix.py index 291cefc7..528028ad 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -113,7 +113,8 @@ def test_build_matrix(ctx_factory, k, curve_fn, op_type, visualize=False): from pytential.symbolic.execution import build_matrix mat = actx.to_numpy( build_matrix(actx, places, sym_op, sym_u, - context=case.knl_concrete_kwargs)) + context=case.knl_concrete_kwargs) + ) if visualize: try: -- GitLab From 9f6f7d92f3a9baae269e8cc0651fbd056dd45550 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 15 Jul 2020 21:26:59 -0500 Subject: [PATCH 04/14] fix comment --- test/extra_matrix_data.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/extra_matrix_data.py b/test/extra_matrix_data.py index bd78f8a9..ea8f238f 100644 --- a/test/extra_matrix_data.py +++ b/test/extra_matrix_data.py @@ -28,8 +28,8 @@ def skeletonization_error(mat, skel, srcindices, p=None): .. math:: \begin{aligned} - \epsilon_{l, ij} = \max_{i \ne j} \|A_{ij} - L_{ii} S_{ij}\|_p, \\ - \epsilon_{r, ij} = \max_{i \ne j} \|A_{ij} - S_{ij} R_{jj}\|_p, \\ + \epsilon_{l, ij} = \|A_{ij} - L_{ii} S_{ij}\|_p, \\ + \epsilon_{r, ij} = \|A_{ij} - S_{ij} R_{jj}\|_p, \\ \epsilon_f = \|A - L S R\|_p. \end{aligned} -- GitLab From 35568a22a1c955d45dc353a3b4ca97cbd894cc99 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 16 Jul 2020 22:06:06 -0500 Subject: [PATCH 05/14] fix some bugs in skeletonization --- pytential/linalg/proxy.py | 20 ++-- pytential/linalg/skeletonization.py | 158 +++++++++++++++------------- pytential/symbolic/matrix.py | 36 +++++-- test/extra_matrix_data.py | 2 +- test/test_linalg_proxy.py | 8 +- 5 files changed, 131 insertions(+), 93 deletions(-) diff --git a/pytential/linalg/proxy.py b/pytential/linalg/proxy.py index 534d8aef..84470c58 100644 --- a/pytential/linalg/proxy.py +++ b/pytential/linalg/proxy.py @@ -62,21 +62,21 @@ def make_block_index(actx, indices, ranges=None): def partition_by_nodes(actx, discr, - tree_kind="adaptive-level-restricted", max_nodes_in_box=None): + tree_kind="adaptive-level-restricted", max_particles_in_box=None): """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 discr: a :class:`meshmode.discretization.Discretization`. :arg tree_kind: if not *None*, it is passed to :class:`boxtree.TreeBuilder`. - :arg max_nodes_in_box: passed to :class:`boxtree.TreeBuilder`. + :arg max_particles_in_box: passed to :class:`boxtree.TreeBuilder`. :return: a :class:`sumpy.tools.BlockIndexRanges`. """ - if max_nodes_in_box is None: + if max_particles_in_box is None: # FIXME: this is just an arbitrary value - max_nodes_in_box = 32 + max_particles_in_box = 32 if tree_kind is not None: from boxtree import box_flags_enum @@ -87,7 +87,7 @@ def partition_by_nodes(actx, discr, from meshmode.dof_array import flatten, thaw tree, _ = builder(actx.queue, flatten(thaw(actx, discr.nodes())), - max_particles_in_box=max_nodes_in_box, + max_particles_in_box=max_particles_in_box, kind=tree_kind) tree = tree.get(actx.queue) @@ -103,7 +103,7 @@ def partition_by_nodes(actx, discr, indices[i] = tree.user_source_ids[box_start:box_end] else: indices = np.arange(0, discr.ndofs, dtype=np.int) - ranges = np.arange(0, discr.ndofs + 1, max_nodes_in_box, dtype=np.int) + ranges = np.arange(0, discr.ndofs + 1, max_particles_in_box, dtype=np.int) return make_block_index(actx, indices, ranges=ranges) @@ -406,7 +406,7 @@ class ProxyGenerator(object): # {{{ gather_block_neighbor_points def gather_block_neighbor_points(actx, discr, indices, pxy, - max_nodes_in_box=None): + max_particles_in_box=None): """Generate a set of neighboring points for each range of points in *discr*. Neighboring points of a range :math:`i` are defined as all the points inside the proxy ball :math:`i` that do not also @@ -419,9 +419,9 @@ def gather_block_neighbor_points(actx, discr, indices, pxy, :return: a :class:`sumpy.tools.BlockIndexRanges`. """ - if max_nodes_in_box is None: + if max_particles_in_box is None: # FIXME: this is a fairly arbitrary value - max_nodes_in_box = 32 + max_particles_in_box = 32 indices = indices.get(actx.queue) @@ -442,7 +442,7 @@ def gather_block_neighbor_points(actx, discr, indices, pxy, from boxtree import TreeBuilder builder = TreeBuilder(actx.context) tree, _ = builder(actx.queue, sources, - max_particles_in_box=max_nodes_in_box) + max_particles_in_box=max_particles_in_box) from boxtree.area_query import AreaQueryBuilder builder = AreaQueryBuilder(actx.context) diff --git a/pytential/linalg/skeletonization.py b/pytential/linalg/skeletonization.py index 95efa15d..c0a756de 100644 --- a/pytential/linalg/skeletonization.py +++ b/pytential/linalg/skeletonization.py @@ -33,6 +33,7 @@ from pytools.obj_array import make_obj_array from pytools import Record from sumpy.tools import MatrixBlockIndexRanges, BlockIndexRanges +from pytential.symbolic.mappers import IdentityMapper, LocationTagger # {{{ helpers @@ -96,22 +97,17 @@ class QBXForcedLimitReplacer(IdentityMapper): return expr.copy(qbx_forced_limit=self.qbx_forced_limit) -class DOFDescReplacer(ToTargetTagger): +class LocationReplacer(LocationTagger): def _default_dofdesc(self, dofdesc): return self.default_where -def prepare_expr(places, expr, auto_where=None, qbx_forced_limit=None): - from pytential.symbolic.execution import _prepare_auto_where - auto_where = _prepare_auto_where(auto_where, places=places) - - from pytential.symbolic.execution import _prepare_expr - expr = _prepare_expr(places, expr, auto_where=auto_where) - - expr = QBXForcedLimitReplacer(qbx_forced_limit)(expr) - expr = DOFDescReplacer(auto_where[0], auto_where[1])(expr) - - return expr +class DOFDescriptorReplacer(LocationReplacer): + def __init__(self, default_source, default_target): + super(DOFDescriptorReplacer, self).__init__( + default_target, default_source=default_source) + self.operand_rec = LocationReplacer( + default_source, default_source=default_source) class BlockEvaluationWrangler: @@ -179,6 +175,15 @@ class BlockEvaluationWrangler: from pytential.symbolic.matrix import FarFieldBlockBuilder self.farfield_block_builder = FarFieldBlockBuilder + def _prepare_farfield_expr(self, places, expr, auto_where=None): + from pytential.symbolic.execution import _prepare_auto_where + auto_where = _prepare_auto_where(auto_where, places=places) + + expr = QBXForcedLimitReplacer(None)(expr) + expr = DOFDescriptorReplacer(auto_where[0], auto_where[1])(expr) + + return expr + def _evaluate(self, actx, places, builder_cls, expr, idomain, index_set, auto_where, **kwargs): domain = self.domains[idomain] @@ -201,7 +206,9 @@ class BlockEvaluationWrangler: def evaluate_source_farfield(self, actx, places, ibrow, ibcol, index_set, auto_where=None): - expr = prepare_expr(places, expr, auto_where=auto_where) + expr = self._prepare_farfield_expr( + places, self.exprs[ibrow], auto_where=auto_where) + return self._evaluate(actx, places, self.farfield_block_builder, expr, ibcol, index_set, auto_where, @@ -210,7 +217,9 @@ class BlockEvaluationWrangler: def evaluate_target_farfield(self, actx, places, ibrow, ibcol, index_set, auto_where=None): - expr = prepare_expr(places, expr, auto_where=auto_where) + expr = self._prepare_farfield_expr( + places, self.exprs[ibrow], auto_where=auto_where) + return self._evaluate(actx, places, self.farfield_block_builder, expr, ibcol, index_set, auto_where, @@ -224,13 +233,19 @@ class BlockEvaluationWrangler: self.exprs[ibrow], ibcol, index_set, auto_where) + def make_block_evaluation_wrangler(places, exprs, input_exprs, - domains=None, context=None, + domains=None, context=None, auto_where=None, _weighted_farfield=None, _farfield_block_builder=None, _nearfield_block_builder=None): - if not is_obj_array(exprs): + from pytential.symbolic.execution import _prepare_auto_where + auto_where = _prepare_auto_where(auto_where, places) + from pytential.symbolic.execution import _prepare_expr + exprs = _prepare_expr(places, exprs, auto_where=auto_where) + + if not (isinstance(exprs, np.ndarray) and exprs.dtype.char == "O"): exprs = make_obj_array([exprs]) try: @@ -238,8 +253,6 @@ def make_block_evaluation_wrangler(places, exprs, input_exprs, except TypeError: input_exprs = [input_exprs] - from pytential.symbolic.execution import _prepare_auto_where - auto_where = _prepare_auto_where(auto_where, places) from pytential.symbolic.execution import _prepare_domains domains = _prepare_domains(len(input_exprs), places, domains, auto_where[0]) @@ -262,11 +275,12 @@ def make_block_evaluation_wrangler(places, exprs, input_exprs, @contextmanager def add_to_geometry_collection(places, proxy): - # NOTE: this is a bit of a hack to keep all the caches in `places` and + # NOTE: this is a giant hack to keep all the caches in `places` and # just add the proxy points to it, since otherwise all the DISCRETIZATION # scope stuff would be recomputed all the time try: + # NOTE: this needs to match `EvaluationMapper.map_common_subexpression` previous_cse_cache = places._get_cache("cse") places.places["proxy"] = proxy yield places @@ -274,14 +288,15 @@ def add_to_geometry_collection(places, proxy): del places.places["proxy"] # NOTE: this is meant to make sure that proxy-related things don't # get cached over multiple runs and lead to some explosion of some sort - places.cache["cse"] = previous_cse_cache + places.caches["cse"] = previous_cse_cache -def make_block_proxy_skeleton(actx, places, proxy, wrangler, indices, +def make_block_proxy_skeleton(actx, places, proxy_generator, wrangler, indices, ibrow, ibcol, source_or_target, max_particles_in_box=None): - """Builds a block matrix that can be used to skeletonize the rows - (targets) of the symbolic matrix block described by ``(ibrow, ibcol)``. + """Builds a block matrix that can be used to skeletonize the + rows (targets) or columns (sources) of the symbolic matrix block + described by ``(ibrow, ibcol)``. """ if source_or_target == "source": @@ -289,11 +304,13 @@ def make_block_proxy_skeleton(actx, places, proxy, wrangler, indices, swap_arg_order = lambda x, y: (x, y) evaluate_farfield = wrangler.evaluate_source_farfield + block_stack = np.vstack elif source_or_target == "target": from pytential.source import PointPotentialSource as ProxyPoints swap_arg_order = lambda x, y: (y, x) evaluate_farfield = wrangler.evaluate_target_farfield + block_stack = np.hstack else: raise ValueError(f"unknown value: '{source_or_target}'") @@ -303,26 +320,28 @@ def make_block_proxy_skeleton(actx, places, proxy, wrangler, indices, dep_source = places.get_geometry(domain.geometry) dep_discr = places.get_discretization(domain.geometry, domain.discr_stage) - pxy = proxy(actx, domain, indices) + pxy = proxy_generator(actx, domain, indices) # }}} # {{{ evaluate (farfield) proxy interactions - with add_to_geometry_collection(places, ProxyPoints(pxy.points)): + with add_to_geometry_collection(places, ProxyPoints(pxy.points)) as pxyplaces: pxyindices = MatrixBlockIndexRanges(actx.context, *swap_arg_order(pxy.indices, indices)) pxymat = evaluate_farfield(actx, pxyplaces, ibrow, ibcol, pxyindices, auto_where=swap_arg_order(domain, "proxy")) if indices.nblocks == 1: - # TODO: evaluate nearfield at the root level? + # NOTE: this doesn't really happen, because this function only gets + # called by skeletonize_block_by_proxy, which already takes care of it return make_block_diag(pxymat, pxyindices.get(actx.queue)) # }}} # {{{ evaluate (nearfield) neighbor interactions + from pytential.linalg.proxy import gather_block_neighbor_points nbrindices = gather_block_neighbor_points( actx, dep_discr, indices, pxy, max_particles_in_box=max_particles_in_box) @@ -335,13 +354,13 @@ def make_block_proxy_skeleton(actx, places, proxy, wrangler, indices, # {{{ concatenate everything to get the blocks ready for ID-ing - pxyindices = pxyindices.get(queue) - nbrindices = nbrindices.get(queue) + pxyindices = pxyindices.get(actx.queue) + nbrindices = nbrindices.get(actx.queue) pxyblk = np.full((indices.nblocks, indices.nblocks), 0, dtype=np.object) for i in range(indices.nblocks): - pxyblk[i, i] = np.hstack(swap_arg_order( - pxyindices.block_take(pxymat, i) + pxyblk[i, i] = block_stack(swap_arg_order( + pxyindices.block_take(pxymat, i), nbrindices.block_take(nbrmat, i), )) # }}} @@ -350,23 +369,9 @@ def make_block_proxy_skeleton(actx, places, proxy, wrangler, indices, def skeletonize_block_by_proxy(actx, - ibcol, ibrow, places, proxy, wrangler, blkindices, + ibcol, ibrow, places, proxy_generator, wrangler, blkindices, id_eps=None, id_rank=None, - tree_max_particles_in_box=None): - r""" - :arg places: a :class:`~meshmode.array_context.ArrayContext`. - :arg proxy: a :class:`~pytential.linalg.proxy.ProxyGenerator`. - :arg wrangler: a :class:`BlockEvaluationWrangler`. - :arg blkindices: a :class:`~sumpy.tools.MatrixBlockIndexRanges`. - - :returns: a tuple ``(L, R, sklindices)`` encoding the block-by-block - decompression of the matrix represented *wrangler*. :math:`L` and - :math:`R` are :math:`n \times n` diagonal block matrix, where - :math:`n` is :attr:`~sumpy.tools.MatrixBlockIndexRanges.nblocks``. The - ``sklindices`` array contains the remaining (skeleton) nodes from - ``blkindices`` after compression. - """ - + max_particles_in_box=None): L = np.full((blkindices.nblocks, blkindices.nblocks), 0, dtype=np.object) R = np.full((blkindices.nblocks, blkindices.nblocks), 0, dtype=np.object) @@ -377,37 +382,34 @@ def skeletonize_block_by_proxy(actx, return L, R, blkindices # construct proxy matrices to skeletonize - src_mat = build_source_skeleton_matrix(queue, - places, proxy, wrangler, blkindices.col, 0, 0, - max_particles_in_box=tree_max_particles_in_box) - tgt_mat = build_target_skeleton_matrix(queue, - places, proxy, wrangler, blkindices.row, 0, 0, - max_particles_in_box=tree_max_particles_in_box) + src_mat = make_block_proxy_skeleton(actx, + places, proxy_generator, wrangler, blkindices.col, + ibrow, ibcol, "source", + max_particles_in_box=max_particles_in_box) + tgt_mat = make_block_proxy_skeleton(actx, + places, proxy_generator, wrangler, blkindices.row, + ibrow, ibcol, "target", + max_particles_in_box=max_particles_in_box) src_skl_indices = np.empty(blkindices.nblocks, dtype=np.object) tgt_skl_indices = np.empty(blkindices.nblocks, dtype=np.object) skl_ranges = np.zeros(blkindices.nblocks + 1, dtype=np.int) - src_indices = blkindices.col.get(queue) - tgt_indices = blkindices.row.get(queue) + src_indices = blkindices.col.get(actx.queue) + tgt_indices = blkindices.row.get(actx.queue) for i in range(blkindices.nblocks): k = id_rank - assert not np.any(np.isnan(src_mat[i, i])), "block {}".format(i) - assert not np.any(np.isinf(src_mat[i, i])), "block {}".format(i) - assert not np.any(np.isnan(tgt_mat[i, i])), "block {}".format(i) - assert not np.any(np.isinf(tgt_mat[i, i])), "block {}".format(i) - # skeletonize target points - k, idx, interp = _interp_decomp(tgt_mat[i, i].T, k, id_eps) + k, idx, interp = interp_decomp(tgt_mat[i, i].T, k, id_eps) assert k > 0 L[i, i] = interp.T tgt_skl_indices[i] = tgt_indices.block_indices(i)[idx[:k]] # skeletonize source points - k, idx, interp = _interp_decomp(src_mat[i, i], k, id_eps) + k, idx, interp = interp_decomp(src_mat[i, i], k, id_eps) assert k > 0 R[i, i] = interp @@ -417,27 +419,39 @@ def skeletonize_block_by_proxy(actx, assert R[i, i].shape == (k, src_mat[i, i].shape[1]) assert L[i, i].shape == (tgt_mat[i, i].shape[0], k) - src_skl_indices = _to_block_index(queue, src_skl_indices, skl_ranges) - tgt_skl_indices = _to_block_index(queue, tgt_skl_indices, skl_ranges) - skl_indices = MatrixBlockIndexRanges(queue.context, - tgt_skl_indices, - src_skl_indices) + from pytential.linalg.proxy import make_block_index + src_skl_indices = make_block_index(actx, np.hstack(src_skl_indices), skl_ranges) + tgt_skl_indices = make_block_index(actx, np.hstack(tgt_skl_indices), skl_ranges) + skl_indices = MatrixBlockIndexRanges(actx.context, + tgt_skl_indices, src_skl_indices) - return L, R, skl_indices + return SkeletonizedBlock(L=L, R=R, sklindices=skl_indices) -def skeletonize_by_proxy(actx, places, proxy, wrangler, blkindices, - id_eps=None, id_rank=None, - max_particles_in_box=None): +def skeletonize_by_proxy(actx, places, proxy_generator, wrangler, blkindices, + id_eps=None, id_rank=None, max_particles_in_box=None): + r""" + :arg places: a :class:`~meshmode.array_context.ArrayContext`. + :arg proxy_generator: a :class:`~pytential.linalg.proxy.ProxyGenerator`. + :arg wrangler: a :class:`BlockEvaluationWrangler`. + :arg blkindices: a :class:`~sumpy.tools.MatrixBlockIndexRanges`. + + :returns: a tuple ``(L, R, sklindices)`` encoding the block-by-block + decompression of the matrix represented *wrangler*. :math:`L` and + :math:`R` are :math:`n \times n` diagonal block matrix, where + :math:`n` is :attr:`~sumpy.tools.MatrixBlockIndexRanges.nblocks``. The + ``sklindices`` array contains the remaining (skeleton) nodes from + ``blkindices`` after compression. + """ nrows = len(wrangler.exprs) ncols = len(wrangler.input_exprs) skel = np.empty((nrows, ncols), dtype=np.object) for ibrow in range(nrows): - for ibcol in range(ncols) + for ibcol in range(ncols): skel[ibrow, ibcol] = skeletonize_block_by_proxy(actx, - ibrow, ibcol, proxy, wrangler, blkindices, + ibrow, ibcol, places, proxy_generator, wrangler, blkindices, id_eps=id_eps, id_rank=id_rank, max_particles_in_box=max_particles_in_box) diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index dc0eb373..70699ce7 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -414,11 +414,13 @@ class MatrixBuilder(MatrixBuilderBase): class P2PMatrixBuilder(MatrixBuilderBase): def __init__(self, actx, dep_expr, other_dep_exprs, - dep_source, dep_discr, places, context, exclude_self=True): + dep_source, dep_discr, places, context, + weightd=False, exclude_self=False): super(P2PMatrixBuilder, self).__init__(actx, dep_expr, other_dep_exprs, dep_source, dep_discr, places, context) + self.weighted = weighted self.exclude_self = exclude_self def map_int_g(self, expr): @@ -459,8 +461,17 @@ class P2PMatrixBuilder(MatrixBuilderBase): targets=flatten(thaw(actx, target_discr.nodes())), sources=flatten(thaw(actx, source_discr.nodes())), **kernel_args) + mat = actx.to_numpy(mat) + + if self.weighted: + from pytential import bind, sym + waa = bind(self.places, sym.weights_and_area_elements( + source_discr.ambient_dim, + dofdesc=expr.source))(self.array_context) + + mat[:, :] *= actx.to_numpy(flatten(waa)) - return actx.to_numpy(mat).dot(rec_density) + return mat.dot(rec_density) # }}} @@ -530,18 +541,22 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase): waa = bind(self.places, sym.weights_and_area_elements( source_discr.ambient_dim, dofdesc=expr.source))(actx) - waa = flatten(waa) + waa = flatten(waa) mat *= waa[self.index_set.linear_col_indices] - return rec_density * actx.to_numpy(mat) + + return actx.to_numpy(mat) * rec_density class FarFieldBlockBuilder(MatrixBlockBuilderBase): def __init__(self, queue, dep_expr, other_dep_exprs, dep_source, dep_discr, - places, index_set, context, exclude_self=False): + places, index_set, context, + weighted=False, exclude_self=False): super(FarFieldBlockBuilder, self).__init__(queue, dep_expr, other_dep_exprs, dep_source, dep_discr, places, index_set, context) + + self.weighted = weighted self.exclude_self = exclude_self def get_dep_variable(self): @@ -593,7 +608,16 @@ class FarFieldBlockBuilder(MatrixBlockBuilderBase): index_set=self.index_set, **kernel_args) - return rec_density * actx.to_numpy(mat) + if self.weighted: + from pytential import bind, sym + waa = bind(self.places, sym.weights_and_area_elements( + source_discr.ambient_dim, + dofdesc=expr.source))(actx) + + waa = flatten(waa) + mat *= waa[self.index_set.linear_col_indices] + + return actx.to_numpy(mat) * rec_density # }}} diff --git a/test/extra_matrix_data.py b/test/extra_matrix_data.py index ea8f238f..163f261b 100644 --- a/test/extra_matrix_data.py +++ b/test/extra_matrix_data.py @@ -130,7 +130,7 @@ class MatrixTestCaseMixin: from pytential.linalg.proxy import partition_by_nodes indices = partition_by_nodes(actx, discr, tree_kind=self.tree_kind, - max_nodes_in_box=max_particles_in_box) + max_particles_in_box=max_particles_in_box) if abs(self.index_sparsity_factor - 1.0) < 1.0e-14: if not matrix_indices: diff --git a/test/test_linalg_proxy.py b/test/test_linalg_proxy.py index fa460141..5aa5f185 100644 --- a/test/test_linalg_proxy.py +++ b/test/test_linalg_proxy.py @@ -388,13 +388,13 @@ def test_skeletonize_by_proxy(ctx_factory, case, visualize=False): queue = cl.CommandQueue(ctx) actx = PyOpenCLArrayContext(queue) - case = case.copy(nblocks=6, op_type="scalar") + case = case.copy(approx_block_count=6, op_type="scalar") logger.info("\n%s", case) # {{{ geometry dd = sym.DOFDescriptor(case.name, discr_stage=case.skel_discr_stage) - qbx = case.get_layer_potential(ctx, case.resolutions[0], case.target_order) + qbx = case.get_layer_potential(actx, case.resolutions[0], case.target_order) places = GeometryCollection(qbx, auto_where=dd) density_discr = places.get_discretization(dd.geometry, dd.discr_stage) @@ -432,7 +432,7 @@ def test_skeletonize_by_proxy(ctx_factory, case, visualize=False): ) # skeleton - from pytential.linalg.skeletonize import skeletonize_by_proxy + from pytential.linalg.skeletonization import skeletonize_by_proxy skeleton = skeletonize_by_proxy(actx, places, proxy, wrangler, srcindices, id_eps=case.id_eps, @@ -444,7 +444,7 @@ def test_skeletonize_by_proxy(ctx_factory, case, visualize=False): # {{{ visualize - if not visualize + if not visualize: return def plot_skeleton(isrc, iskl, name): -- GitLab From 7819ce71b1f627b6215e6562abdc48d14f472793 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 22 Jul 2020 15:59:34 -0500 Subject: [PATCH 06/14] fix skeletonization The LocationTagger wasn't tagging IntGs properly because they don't use _default_dofdesc in the base class. Overwriting map_int_g to correctly replace the source and target seems to have worked --- pytential/linalg/skeletonization.py | 89 ++++++++++++++++++++--------- pytential/symbolic/matrix.py | 10 ++-- pytential/utils.py | 2 +- test/extra_matrix_data.py | 12 ++-- test/test_linalg_proxy.py | 56 ++++++++++++++---- 5 files changed, 119 insertions(+), 50 deletions(-) diff --git a/pytential/linalg/skeletonization.py b/pytential/linalg/skeletonization.py index c0a756de..9aa105ea 100644 --- a/pytential/linalg/skeletonization.py +++ b/pytential/linalg/skeletonization.py @@ -76,13 +76,6 @@ def make_block_diag(blk, blkindices): return diag -class SkeletonizedBlock(Record): - """ - .. attribute:: L - .. attribute:: R - .. attribute:: sklindices - """ - # }}} @@ -101,6 +94,17 @@ class LocationReplacer(LocationTagger): def _default_dofdesc(self, dofdesc): return self.default_where + def map_int_g(self, expr): + return type(expr)( + expr.kernel, + self.operand_rec(expr.density), + expr.qbx_forced_limit, + self.default_source, self.default_where, + kernel_arguments=dict( + (name, self.operand_rec(arg_expr)) + for name, arg_expr in expr.kernel_arguments.items() + )) + class DOFDescriptorReplacer(LocationReplacer): def __init__(self, default_source, default_target): @@ -179,7 +183,7 @@ class BlockEvaluationWrangler: from pytential.symbolic.execution import _prepare_auto_where auto_where = _prepare_auto_where(auto_where, places=places) - expr = QBXForcedLimitReplacer(None)(expr) + expr = QBXForcedLimitReplacer(qbx_forced_limit=None)(expr) expr = DOFDescriptorReplacer(auto_where[0], auto_where[1])(expr) return expr @@ -271,7 +275,7 @@ def make_block_evaluation_wrangler(places, exprs, input_exprs, # }}} -# {{{ skeletonize_by_proxy +# {{{ skeletonize_block_by_proxy @contextmanager def add_to_geometry_collection(places, proxy): @@ -280,19 +284,26 @@ def add_to_geometry_collection(places, proxy): # scope stuff would be recomputed all the time try: - # NOTE: this needs to match `EvaluationMapper.map_common_subexpression` - previous_cse_cache = places._get_cache("cse") - places.places["proxy"] = proxy - yield places + pxyplaces = places.merge({"proxy": proxy}) + yield pxyplaces finally: - del places.places["proxy"] - # NOTE: this is meant to make sure that proxy-related things don't - # get cached over multiple runs and lead to some explosion of some sort - places.caches["cse"] = previous_cse_cache - - -def make_block_proxy_skeleton(actx, places, proxy_generator, wrangler, indices, - ibrow, ibcol, source_or_target, + pass + + # try: + # # NOTE: this needs to match `EvaluationMapper.map_common_subexpression` + # previous_cse_cache = places._get_cache("cse") + # places.places["proxy"] = proxy + # yield places + # finally: + # del places.places["proxy"] + # # NOTE: this is meant to make sure that proxy-related things don't + # # get cached over multiple runs and lead to some explosion of some sort + # places.caches["cse"] = previous_cse_cache + + +def make_block_proxy_skeleton(actx, ibrow, ibcol, + places, proxy_generator, wrangler, indices, + source_or_target=None, max_particles_in_box=None): """Builds a block matrix that can be used to skeletonize the rows (targets) or columns (sources) of the symbolic matrix block @@ -382,13 +393,13 @@ def skeletonize_block_by_proxy(actx, return L, R, blkindices # construct proxy matrices to skeletonize - src_mat = make_block_proxy_skeleton(actx, + src_mat = make_block_proxy_skeleton(actx, ibrow, ibcol, places, proxy_generator, wrangler, blkindices.col, - ibrow, ibcol, "source", + source_or_target="source", max_particles_in_box=max_particles_in_box) - tgt_mat = make_block_proxy_skeleton(actx, + tgt_mat = make_block_proxy_skeleton(actx, ibrow, ibcol, places, proxy_generator, wrangler, blkindices.row, - ibrow, ibcol, "target", + source_or_target="target", max_particles_in_box=max_particles_in_box) src_skl_indices = np.empty(blkindices.nblocks, dtype=np.object) @@ -402,6 +413,11 @@ def skeletonize_block_by_proxy(actx, k = id_rank # skeletonize target points + assert not np.any(np.isinf(tgt_mat[i, i])), \ + np.where(np.isinf(tgt_mat[i, i])) + assert not np.any(np.isnan(tgt_mat[i, i])), \ + np.where(np.isnan(tgt_mat[i, i])) + k, idx, interp = interp_decomp(tgt_mat[i, i].T, k, id_eps) assert k > 0 @@ -409,6 +425,11 @@ def skeletonize_block_by_proxy(actx, tgt_skl_indices[i] = tgt_indices.block_indices(i)[idx[:k]] # skeletonize source points + assert not np.any(np.isinf(src_mat[i, i])), \ + np.where(np.isinf(src_mat[i, i])) + assert not np.any(np.isnan(src_mat[i, i])), \ + np.where(np.isnan(src_mat[i, i])) + k, idx, interp = interp_decomp(src_mat[i, i], k, id_eps) assert k > 0 @@ -425,7 +446,19 @@ def skeletonize_block_by_proxy(actx, skl_indices = MatrixBlockIndexRanges(actx.context, tgt_skl_indices, src_skl_indices) - return SkeletonizedBlock(L=L, R=R, sklindices=skl_indices) + return L, R, skl_indices + +# }}} + + +# {{{ skeletonize_by_proxy + +class SkeletonizedBlock(Record): + """ + .. attribute:: L + .. attribute:: R + .. attribute:: sklindices + """ def skeletonize_by_proxy(actx, places, proxy_generator, wrangler, blkindices, @@ -450,11 +483,13 @@ def skeletonize_by_proxy(actx, places, proxy_generator, wrangler, blkindices, for ibrow in range(nrows): for ibcol in range(ncols): - skel[ibrow, ibcol] = skeletonize_block_by_proxy(actx, + L, R, sklindices = skeletonize_block_by_proxy(actx, ibrow, ibcol, places, proxy_generator, wrangler, blkindices, id_eps=id_eps, id_rank=id_rank, max_particles_in_box=max_particles_in_box) + skel[ibrow, ibcol] = SkeletonizedBlock(L=L, R=R, sklindices=sklindices) + return skel # }}} diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 70699ce7..f4457e06 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -572,9 +572,6 @@ class FarFieldBlockBuilder(MatrixBlockBuilderBase): target_discr = self.places.get_discretization( expr.target.geometry, expr.target.discr_stage) - if source_discr is not target_discr: - raise NotImplementedError - rec_density = self._blk_mapper.rec(expr.density) if is_zero(rec_density): return 0 @@ -601,10 +598,10 @@ class FarFieldBlockBuilder(MatrixBlockBuilderBase): mat_gen = P2PMatrixBlockGenerator(actx.context, (kernel,), exclude_self=self.exclude_self) - from meshmode.dof_array import flatten, thaw + from pytential.utils import flatten_if_needed _, (mat,) = mat_gen(actx.queue, - targets=flatten(thaw(actx, target_discr.nodes())), - sources=flatten(thaw(actx, source_discr.nodes())), + targets=flatten_if_needed(actx, target_discr.nodes()), + sources=flatten_if_needed(actx, source_discr.nodes()), index_set=self.index_set, **kernel_args) @@ -614,6 +611,7 @@ class FarFieldBlockBuilder(MatrixBlockBuilderBase): source_discr.ambient_dim, dofdesc=expr.source))(actx) + from meshmode.dof_array import flatten waa = flatten(waa) mat *= waa[self.index_set.linear_col_indices] diff --git a/pytential/utils.py b/pytential/utils.py index f1e9f0d1..88fb3ad3 100644 --- a/pytential/utils.py +++ b/pytential/utils.py @@ -36,7 +36,7 @@ def flatten_if_needed(actx: PyOpenCLArrayContext, ary: np.ndarray): return obj_array_vectorize_n_args(flatten_if_needed, actx, ary) if not isinstance(ary, DOFArray): - return ary + return ary.copy(actx.queue) if ary.array_context is None: ary = thaw(actx, ary) diff --git a/test/extra_matrix_data.py b/test/extra_matrix_data.py index 163f261b..805bc3ef 100644 --- a/test/extra_matrix_data.py +++ b/test/extra_matrix_data.py @@ -58,8 +58,8 @@ def skeletonization_error(mat, skel, srcindices, p=None): return blk.indices[blk.ranges[i]:blk.ranges[i + 1]] # compute max block error - err_l = np.empty((nblocks, nblocks)) - err_l = np.empty((nblocks, nblocks)) + err_l = np.zeros((nblocks, nblocks)) + err_r = np.zeros((nblocks, nblocks)) for i, j in product(range(nblocks), repeat=2): if i == j: continue @@ -82,10 +82,10 @@ def skeletonization_error(mat, skel, srcindices, p=None): # compute full matrix error from pytential.symbolic.execution import _bmat - A = _bmat(A, dtype=mat.dtype) - L = _bmat(L, dtype=mat.dtype) - S = _bmat(S, dtype=mat.dtype) - R = _bmat(R, dtype=mat.dtype) + A = _bmat(A, dtypes=[mat.dtype]) + L = _bmat(L, dtypes=[mat.dtype]) + S = _bmat(S, dtypes=[mat.dtype]) + R = _bmat(R, dtypes=[mat.dtype]) assert L.shape == (A.shape[0], S.shape[0]) assert R.shape == (S.shape[1], A.shape[1]) diff --git a/test/test_linalg_proxy.py b/test/test_linalg_proxy.py index 5aa5f185..eac94e01 100644 --- a/test/test_linalg_proxy.py +++ b/test/test_linalg_proxy.py @@ -31,9 +31,10 @@ import pyopencl as cl from pytential import bind, sym from pytential import GeometryCollection +from pytential.utils import flatten_to_numpy from meshmode.array_context import PyOpenCLArrayContext -from meshmode.mesh.generation import ellipse +from meshmode.mesh.generation import ellipse, NArmedStarfish import pytest from pyopencl.tools import ( # noqa @@ -66,7 +67,6 @@ def plot_partition_indices(actx, discr, indices, **kwargs): pt.savefig("test_partition_{1}_{3}d_ranges_{2}.png".format(*args)) pt.clf() - from pytential.utils import flatten_to_numpy if discr.ambient_dim == 2: sources = flatten_to_numpy(actx, discr.nodes()) @@ -108,6 +108,11 @@ PROXY_TEST_CASES = [ name="ellipse", target_order=7, curve_fn=partial(ellipse, 3.0)), + extra.CurveTestCase( + name="starfish", + target_order=4, + curve_fn=NArmedStarfish(5, 0.25), + resolutions=[32]), extra.TorusTestCase( target_order=2, resolutions=[0]) @@ -172,7 +177,6 @@ def test_proxy_generator(ctx_factory, case, index_sparsity_factor, visualize=Fal from pytential.linalg.proxy import ProxyGenerator proxies = ProxyGenerator(places)(actx, places.auto_source, srcindices) - from pytential.utils import flatten_to_numpy proxies = proxies.get(queue) pxypoints = np.vstack(proxies.points) pxycenters = np.vstack(proxies.centers) @@ -305,7 +309,6 @@ def test_neighbor_points(ctx_factory, case, index_sparsity_factor, visualize=Fal srcindices = srcindices.get(queue) nbrindices = nbrindices.get(queue) - from pytential.utils import flatten_to_numpy proxies = proxies.get(queue) pxypoints = np.vstack(proxies.points) pxycenters = np.vstack(proxies.centers) @@ -408,7 +411,7 @@ def test_skeletonize_by_proxy(ctx_factory, case, visualize=False): from pytential.linalg.proxy import ProxyGenerator from pytential.linalg.skeletonization import make_block_evaluation_wrangler - proxy = ProxyGenerator(places, + proxy_generator = ProxyGenerator(places, radius_factor=case.proxy_radius_factor, approx_nproxy=case.proxy_approx_count) @@ -432,13 +435,24 @@ def test_skeletonize_by_proxy(ctx_factory, case, visualize=False): ) # skeleton - from pytential.linalg.skeletonization import skeletonize_by_proxy - skeleton = skeletonize_by_proxy(actx, - places, proxy, wrangler, srcindices, + from pytential.linalg.skeletonization import skeletonize_block_by_proxy + L, R, sklindices = skeletonize_block_by_proxy(actx, 0, 0, + places, proxy_generator, wrangler, srcindices, id_eps=case.id_eps, max_particles_in_box=case.max_particles_in_box) - err_l, err_r, err_f = extra.skeletonization_error(mat, skeleton, srcindices) + from pytential.linalg.skeletonization import SkeletonizedBlock + skeleton = SkeletonizedBlock(L=L, R=R, sklindices=sklindices.get(queue)) + + err_l, err_r, err_f = extra.skeletonization_error( + mat, skeleton, srcindices.get(queue)) + + logger.info("error: L %.5e R %.5e F %.5e", + la.norm(err_l, np.inf), la.norm(err_r, np.inf), err_f) + + assert la.norm(err_l, np.inf) < 5.0e+1 * case.id_eps + assert la.norm(err_r, np.inf) < 5.0e+1 * case.id_eps + assert err_f < 5.0e+1 * case.id_eps # }}} @@ -447,18 +461,40 @@ def test_skeletonize_by_proxy(ctx_factory, case, visualize=False): if not visualize: return + import matplotlib.pyplot as pt + pt.imshow(np.log10(err_l + 1.0e-16)) + pt.colorbar() + pt.savefig("test_skeletonize_by_proxy_err_l") + pt.clf() + + pt.imshow(np.log10(err_r + 1.0e-16)) + pt.colorbar() + pt.savefig("test_skeletonize_by_proxy_err_r") + pt.clf() + def plot_skeleton(isrc, iskl, name): pt.figure(figsize=(10, 10), dpi=300) pt.plot(sources[0][isrc.indices], sources[1][isrc.indices], 'ko', alpha=0.5) - for i in range(blkindices.nblocks): + ax = pt.gca() + for i in range(srcindices.nblocks): iblk = iskl.block_indices(i) pt.plot(sources[0][iblk], sources[1][iblk], 'o') + c = pt.Circle(pxycenters[:, i], pxyradii[i], color='k', alpha=0.1) + ax.add_artist(c) + + ax.set_aspect("equal") + pt.xlim([-1.5, 1.5]) + pt.ylim([-1.5, 1.5]) pt.savefig(f"test_skeletonize_by_proxy_{name}") pt.clf() + pxy = proxy_generator(actx, wrangler.domains[0], srcindices.row).get(queue) + pxycenters = np.vstack(pxy.centers) + pxyradii = pxy.radii + if places.ambient_dim == 2: sources = flatten_to_numpy(actx, density_discr.nodes()) srcindices = srcindices.get(queue) -- GitLab From a6ba288f96101068e9d5d50d68f0c5c7def6daa7 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 23 Jul 2020 09:48:55 -0500 Subject: [PATCH 07/14] flake8 fixes --- pytential/linalg/proxy.py | 2 +- pytential/linalg/skeletonization.py | 14 ++++---------- pytential/symbolic/matrix.py | 2 +- test/test_linalg_proxy.py | 1 - 4 files changed, 6 insertions(+), 13 deletions(-) diff --git a/pytential/linalg/proxy.py b/pytential/linalg/proxy.py index 84470c58..d90cfa18 100644 --- a/pytential/linalg/proxy.py +++ b/pytential/linalg/proxy.py @@ -27,7 +27,7 @@ import numpy as np import numpy.linalg as la from pytools.obj_array import make_obj_array -from pytools import memoize_method, memoize_in +from pytools import memoize_method from sumpy.tools import BlockIndexRanges from boxtree.tools import DeviceDataRecord diff --git a/pytential/linalg/skeletonization.py b/pytential/linalg/skeletonization.py index 9aa105ea..e2d065db 100644 --- a/pytential/linalg/skeletonization.py +++ b/pytential/linalg/skeletonization.py @@ -25,14 +25,11 @@ THE SOFTWARE. from contextlib import contextmanager import numpy as np -import numpy.linalg as la - -import pyopencl as cl from pytools.obj_array import make_obj_array from pytools import Record -from sumpy.tools import MatrixBlockIndexRanges, BlockIndexRanges +from sumpy.tools import MatrixBlockIndexRanges from pytential.symbolic.mappers import IdentityMapper, LocationTagger @@ -126,7 +123,7 @@ class BlockEvaluationWrangler: weighted_farfield=None, farfield_block_builder=None, nearfield_block_builder=None): - """ + r""" :arg exprs: an :class:`~numpy.ndarray` of expressions (layer potentials) that correspond to the output blocks of the matrix. :arg input_exprs: a :class:`list` of densities that correspond to the @@ -237,7 +234,6 @@ class BlockEvaluationWrangler: self.exprs[ibrow], ibcol, index_set, auto_where) - def make_block_evaluation_wrangler(places, exprs, input_exprs, domains=None, context=None, auto_where=None, _weighted_farfield=None, @@ -313,13 +309,13 @@ def make_block_proxy_skeleton(actx, ibrow, ibcol, if source_or_target == "source": from pytential.target import PointsTarget as ProxyPoints - swap_arg_order = lambda x, y: (x, y) + swap_arg_order = lambda x, y: (x, y) # noqa: E731 evaluate_farfield = wrangler.evaluate_source_farfield block_stack = np.vstack elif source_or_target == "target": from pytential.source import PointPotentialSource as ProxyPoints - swap_arg_order = lambda x, y: (y, x) + swap_arg_order = lambda x, y: (y, x) # noqa: E731 evaluate_farfield = wrangler.evaluate_target_farfield block_stack = np.hstack else: @@ -328,9 +324,7 @@ def make_block_proxy_skeleton(actx, ibrow, ibcol, # {{{ generate proxies domain = wrangler.domains[ibcol] - dep_source = places.get_geometry(domain.geometry) dep_discr = places.get_discretization(domain.geometry, domain.discr_stage) - pxy = proxy_generator(actx, domain, indices) # }}} diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index f4457e06..b6ccb5de 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -415,7 +415,7 @@ class MatrixBuilder(MatrixBuilderBase): class P2PMatrixBuilder(MatrixBuilderBase): def __init__(self, actx, dep_expr, other_dep_exprs, dep_source, dep_discr, places, context, - weightd=False, exclude_self=False): + weighted=False, exclude_self=False): super(P2PMatrixBuilder, self).__init__(actx, dep_expr, other_dep_exprs, dep_source, dep_discr, places, context) diff --git a/test/test_linalg_proxy.py b/test/test_linalg_proxy.py index eac94e01..839b4074 100644 --- a/test/test_linalg_proxy.py +++ b/test/test_linalg_proxy.py @@ -310,7 +310,6 @@ def test_neighbor_points(ctx_factory, case, index_sparsity_factor, visualize=Fal nbrindices = nbrindices.get(queue) proxies = proxies.get(queue) - pxypoints = np.vstack(proxies.points) pxycenters = np.vstack(proxies.centers) nodes = np.vstack(flatten_to_numpy(actx, density_discr.nodes())) -- GitLab From 2f94dfe9a996e663dfed98d9db196c9bdb358072 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 23 Jul 2020 09:50:48 -0500 Subject: [PATCH 08/14] remove copy --- pytential/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pytential/utils.py b/pytential/utils.py index 88fb3ad3..f1e9f0d1 100644 --- a/pytential/utils.py +++ b/pytential/utils.py @@ -36,7 +36,7 @@ def flatten_if_needed(actx: PyOpenCLArrayContext, ary: np.ndarray): return obj_array_vectorize_n_args(flatten_if_needed, actx, ary) if not isinstance(ary, DOFArray): - return ary.copy(actx.queue) + return ary if ary.array_context is None: ary = thaw(actx, ary) -- GitLab From 06285a0a203d1f83fafd3674830899b914cb9eea Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 23 Jul 2020 14:15:40 -0500 Subject: [PATCH 09/14] make proxy coordinates indepedent arrays Before they were slices of a bigger array and had an `offset`, which caused problems with several kernels. --- pytential/linalg/proxy.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pytential/linalg/proxy.py b/pytential/linalg/proxy.py index d90cfa18..d6927caa 100644 --- a/pytential/linalg/proxy.py +++ b/pytential/linalg/proxy.py @@ -384,10 +384,10 @@ class ProxyGenerator(object): pxy_nr_base += self.nproxy proxies = make_obj_array([ - actx.freeze(p) for p in actx.from_numpy(proxies) + actx.freeze(actx.from_numpy(p)) for p in proxies ]) centers = make_obj_array([ - actx.freeze(c) for c in centers_dev + actx.freeze(actx.from_numpy(c)) for c in centers ]) pxyindices = np.arange(0, nproxy, dtype=indices.indices.dtype) pxyranges = np.arange(0, nproxy + 1, self.nproxy) -- GitLab From 1746a2d70993acb0acc43f903d4b97c5a762d15c Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 23 Jul 2020 18:27:02 -0500 Subject: [PATCH 10/14] increase assert tolerance in 3d --- test/test_linalg_proxy.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/test/test_linalg_proxy.py b/test/test_linalg_proxy.py index 839b4074..1f07a7e7 100644 --- a/test/test_linalg_proxy.py +++ b/test/test_linalg_proxy.py @@ -449,9 +449,12 @@ def test_skeletonize_by_proxy(ctx_factory, case, visualize=False): logger.info("error: L %.5e R %.5e F %.5e", la.norm(err_l, np.inf), la.norm(err_r, np.inf), err_f) - assert la.norm(err_l, np.inf) < 5.0e+1 * case.id_eps - assert la.norm(err_r, np.inf) < 5.0e+1 * case.id_eps - assert err_f < 5.0e+1 * case.id_eps + # FIXME: why is the 3D error so large? + rtol = 5 * 10**places.ambient_dim * case.id_eps + + assert la.norm(err_l, np.inf) < rtol + assert la.norm(err_r, np.inf) < rtol + assert err_f < rtol # }}} -- GitLab From 3d65d223d0bfc230f167857cd98a9dae8a37983e Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 13 Sep 2020 11:30:07 -0500 Subject: [PATCH 11/14] fix merge --- test/test_linalg_proxy.py | 49 --------------------------------------- 1 file changed, 49 deletions(-) diff --git a/test/test_linalg_proxy.py b/test/test_linalg_proxy.py index 2c109dea..65063752 100644 --- a/test/test_linalg_proxy.py +++ b/test/test_linalg_proxy.py @@ -373,55 +373,6 @@ def test_neighbor_points(ctx_factory, case, index_sparsity_factor, visualize=Fal ("marker", marker_dev), ]) - density_nodes = flatten_to_numpy(actx, density_discr.nodes()) - nodes = flatten_to_numpy(actx, nodes) - ranges = actx.to_numpy(ranges) - - for i in range(srcindices.nblocks): - isrc = srcindices.block_indices(i) - inbr = nbrindices.block_indices(i) - iall = np.s_[ranges[i]:ranges[i + 1]] - - pt.figure(figsize=(10, 8)) - pt.plot(density_nodes[0], density_nodes[1], - "ko", ms=2.0, alpha=0.5) - pt.plot(density_nodes[0][srcindices.indices], - density_nodes[1][srcindices.indices], - "o", ms=2.0) - pt.plot(density_nodes[0][isrc], density_nodes[1][isrc], - "o", ms=2.0) - pt.plot(density_nodes[0][inbr], density_nodes[1][inbr], - "o", ms=2.0) - pt.plot(nodes[0][iall], nodes[1][iall], - "x", ms=2.0) - pt.xlim([-1.5, 1.5]) - pt.ylim([-1.5, 1.5]) - - filename = f"test_area_query_{ambient_dim}d_{i:04}.png" - pt.savefig(filename, dpi=300) - pt.clf() - elif ambient_dim == 3: - from meshmode.discretization.visualization import make_visualizer - marker = np.empty(density_discr.ndofs) - - for i in range(srcindices.nblocks): - isrc = srcindices.block_indices(i) - inbr = nbrindices.block_indices(i) - - marker.fill(0.0) - marker[srcindices.indices] = 0.0 - marker[isrc] = -42.0 - marker[inbr] = +42.0 - - from meshmode.dof_array import unflatten - marker_dev = unflatten(actx, density_discr, actx.from_numpy(marker)) - - vis = make_visualizer(actx, density_discr, 10) - filename = f"test_area_query_{ambient_dim}d_{i:04}.vtu" - vis.write_vtk_file(filename, [ - ("marker", marker_dev), - ]) ->>>>>>> master # }}} # }}} -- GitLab From 54c3dd4fa31fe04e368221327c5b9603c3cce3da Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 13 Sep 2020 12:05:29 -0500 Subject: [PATCH 12/14] linting fixes --- .gitlab-ci.yml | 2 +- pytential/linalg/proxy.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 81b69ec1..308375dc 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -88,7 +88,7 @@ Pylint: # Pylint won't find the Cython bits without this - PROJECT_INSTALL_FLAGS="--editable" - export PY_EXE=python3 - - EXTRA_INSTALL="Cython pybind11 numpy mako matplotlib" + - EXTRA_INSTALL="Cython pybind11 numpy scipy mako matplotlib" - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/prepare-and-run-pylint.sh - ". ./prepare-and-run-pylint.sh pytential test/test_*.py" tags: diff --git a/pytential/linalg/proxy.py b/pytential/linalg/proxy.py index 26fd740a..7450e3d3 100644 --- a/pytential/linalg/proxy.py +++ b/pytential/linalg/proxy.py @@ -39,10 +39,10 @@ Proxy Point Generation ~~~~~~~~~~~~~~~~~~~~~~ .. autoclass:: ProxyGenerator +.. autoclass:: BlockProxyPoints .. autofunction:: partition_by_nodes .. autofunction:: gather_block_neighbor_points - """ -- GitLab From 78130050aadf3d9f1637dfa3b67c90c9fac7f83d Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 16 Sep 2020 14:26:51 -0500 Subject: [PATCH 13/14] tests that the id_decomp reaches tolerance --- pytential/linalg/skeletonization.py | 17 +++++++- test/test_linalg_proxy.py | 65 ++++++++++++++++++++++------- 2 files changed, 66 insertions(+), 16 deletions(-) diff --git a/pytential/linalg/skeletonization.py b/pytential/linalg/skeletonization.py index e2d065db..131aca97 100644 --- a/pytential/linalg/skeletonization.py +++ b/pytential/linalg/skeletonization.py @@ -373,7 +373,7 @@ def make_block_proxy_skeleton(actx, ibrow, ibcol, return pxyblk -def skeletonize_block_by_proxy(actx, +def _skeletonize_block_by_proxy_with_mats(actx, ibcol, ibrow, places, proxy_generator, wrangler, blkindices, id_eps=None, id_rank=None, max_particles_in_box=None): @@ -440,7 +440,20 @@ def skeletonize_block_by_proxy(actx, skl_indices = MatrixBlockIndexRanges(actx.context, tgt_skl_indices, src_skl_indices) - return L, R, skl_indices + return L, R, skl_indices, src_mat, tgt_mat + + +def skeletonize_block_by_proxy(actx, + ibcol, ibrow, places, proxy_generator, wrangler, blkindices, + id_eps=None, id_rank=None, + max_particles_in_box=None): + L, R, sklindices, _, _ = _skeletonize_block_by_proxy_with_mats(actx, + ibcol, ibrow, places, proxy_generator, wrangler, blkindices, + id_eps=id_eps, + id_rank=id_rank, + max_particles_in_box=max_particles_in_box) + + return L, R, sklindices # }}} diff --git a/test/test_linalg_proxy.py b/test/test_linalg_proxy.py index 65063752..629e221e 100644 --- a/test/test_linalg_proxy.py +++ b/test/test_linalg_proxy.py @@ -388,7 +388,7 @@ def test_skeletonize_by_proxy(ctx_factory, case, visualize=False): queue = cl.CommandQueue(ctx) actx = PyOpenCLArrayContext(queue) - case = case.copy(approx_block_count=6, op_type="scalar") + case = case.copy(approx_block_count=6, op_type="scalar", id_eps=1.0e-8) logger.info("\n%s", case) # {{{ geometry @@ -422,7 +422,7 @@ def test_skeletonize_by_proxy(ctx_factory, case, visualize=False): # }}} - # {{{ check skeletonize + # {{{ check proxy id decomposition # dense matrix from pytential.symbolic.execution import build_matrix @@ -432,26 +432,63 @@ def test_skeletonize_by_proxy(ctx_factory, case, visualize=False): ) # skeleton - from pytential.linalg.skeletonization import skeletonize_block_by_proxy - L, R, sklindices = skeletonize_block_by_proxy(actx, 0, 0, - places, proxy_generator, wrangler, srcindices, + from pytential.linalg.skeletonization import \ + _skeletonize_block_by_proxy_with_mats + L, R, sklindices, src, tgt = _skeletonize_block_by_proxy_with_mats(actx, + 0, 0, places, proxy_generator, wrangler, srcindices, id_eps=case.id_eps, max_particles_in_box=case.max_particles_in_box) + srcindices = srcindices.get(queue) + sklindices = sklindices.get(queue) + for i in range(sklindices.nblocks): + # targets (rows) + bi = np.searchsorted( + srcindices.row.block_indices(i), + sklindices.row.block_indices(i), + ) + + A = tgt[i, i] + S = A[bi, :] + tgt_error = la.norm(A - L[i, i] @ S) / la.norm(A) + + # sources (columns) + bj = np.searchsorted( + srcindices.col.block_indices(i), + sklindices.col.block_indices(i), + ) + + A = src[i, i] + S = A[:, bj] + src_error = la.norm(A - S @ R[i, i]) / la.norm(A) + + logger.info("[%04d] id_eps %.5e src %.5e tgt %.5e rank %d/%d", + i, case.id_eps, + src_error, tgt_error, R[i, i].shape[0], R[i, i].shape[1]) + + assert src_error < 6 * case.id_eps + assert tgt_error < 6 * case.id_eps + + # }}} + + # {{{ check skeletonize + from pytential.linalg.skeletonization import SkeletonizedBlock - skeleton = SkeletonizedBlock(L=L, R=R, sklindices=sklindices.get(queue)) + skeleton = SkeletonizedBlock(L=L, R=R, sklindices=sklindices) - err_l, err_r, err_f = extra.skeletonization_error( - mat, skeleton, srcindices.get(queue)) + blk_err_l, blk_err_r, err_f = extra.skeletonization_error( + mat, skeleton, srcindices) + err_l = la.norm(blk_err_l, np.inf) + err_r = la.norm(blk_err_r, np.inf) - logger.info("error: L %.5e R %.5e F %.5e", - la.norm(err_l, np.inf), la.norm(err_r, np.inf), err_f) + logger.info("error: id_eps %.5e L %.5e R %.5e F %.5e", + case.id_eps, err_l, err_r, err_f) # FIXME: why is the 3D error so large? rtol = 5 * 10**places.ambient_dim * case.id_eps - assert la.norm(err_l, np.inf) < rtol - assert la.norm(err_r, np.inf) < rtol + assert err_l < rtol + assert err_r < rtol assert err_f < rtol # }}} @@ -462,12 +499,12 @@ def test_skeletonize_by_proxy(ctx_factory, case, visualize=False): return import matplotlib.pyplot as pt - pt.imshow(np.log10(err_l + 1.0e-16)) + pt.imshow(np.log10(blk_err_l + 1.0e-16)) pt.colorbar() pt.savefig("test_skeletonize_by_proxy_err_l") pt.clf() - pt.imshow(np.log10(err_r + 1.0e-16)) + pt.imshow(np.log10(blk_err_r + 1.0e-16)) pt.colorbar() pt.savefig("test_skeletonize_by_proxy_err_r") pt.clf() -- GitLab From 398771048cf5f45ee070e449e5f68479aad22248 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 17 Sep 2020 14:09:29 -0500 Subject: [PATCH 14/14] set proxy center in middle of block bounding box --- pytential/linalg/proxy.py | 118 +++++++++++++++++++++----------------- test/test_linalg_proxy.py | 88 ++++++++++++++++++---------- 2 files changed, 125 insertions(+), 81 deletions(-) diff --git a/pytential/linalg/proxy.py b/pytential/linalg/proxy.py index 7450e3d3..0526329f 100644 --- a/pytential/linalg/proxy.py +++ b/pytential/linalg/proxy.py @@ -250,44 +250,54 @@ class ProxyGenerator: radius_expr = "{f} * rqbx" radius_expr = radius_expr.format(f=self.radius_factor) - # NOTE: centers of mass are computed using a second-order approximation knl = lp.make_kernel([ "{[irange]: 0 <= irange < nranges}", - "{[i]: 0 <= i < npoints}", - "{[idim]: 0 <= idim < dim}" + "{[i]: 0 <= i < npoints}", "{[j]: 0 <= j < npoints}", + "{[idim]: 0 <= idim < ndim}" ], - [""" + """ for irange <> ioffset = srcranges[irange] <> npoints = srcranges[irange + 1] - srcranges[irange] - proxy_center[idim, irange] = 1.0 / npoints * \ - reduce(sum, i, sources[idim, srcindices[i + ioffset]]) \ - {{dup=idim:i}} + for idim + <> bbox_max = \ + simul_reduce(max, i, sources[idim, srcindices[i + ioffset]]) + <> bbox_min = \ + simul_reduce(min, i, sources[idim, srcindices[i + ioffset]]) - <> rblk = simul_reduce(max, i, sqrt(simul_reduce(sum, idim, \ - (proxy_center[idim, irange] - - sources[idim, srcindices[i + ioffset]]) ** 2))) + proxy_center[idim, irange] = (bbox_max + bbox_min) / 2.0 + end - <> rqbx_int = simul_reduce(max, i, sqrt(simul_reduce(sum, idim, \ + <> rblk = simul_reduce(max, j, sqrt(simul_reduce(sum, idim, \ + (proxy_center[idim, irange] - + sources[idim, srcindices[j + ioffset]]) ** 2))) \ + {dup=idim} + """ + + (""" + <> rqbx_int = simul_reduce(max, j, sqrt(simul_reduce(sum, idim, \ (proxy_center[idim, irange] - - center_int[idim, srcindices[i + ioffset]]) ** 2)) + \ - expansion_radii[srcindices[i + ioffset]]) - <> rqbx_ext = simul_reduce(max, i, sqrt(simul_reduce(sum, idim, \ + center_int[idim, srcindices[j + ioffset]]) ** 2)) + \ + expansion_radii[srcindices[j + ioffset]]) \ + {dup=idim} + <> rqbx_ext = simul_reduce(max, j, sqrt(simul_reduce(sum, idim, \ (proxy_center[idim, irange] - - center_ext[idim, srcindices[i + ioffset]]) ** 2)) + \ - expansion_radii[srcindices[i + ioffset]]) + center_ext[idim, srcindices[j + ioffset]]) ** 2)) + \ + expansion_radii[srcindices[j + ioffset]]) \ + {dup=idim} <> rqbx = rqbx_int if rqbx_ext < rqbx_int else rqbx_ext - + """ if self.radius_factor > 1.0e-14 else "<> rqbx = 0.0") + + """ proxy_radius[irange] = {radius_expr} end - """.format(radius_expr=radius_expr)], - [ - lp.GlobalArg("sources", None, - shape=(self.ambient_dim, "nsources"), dim_tags="sep,C"), + """.format(radius_expr=radius_expr), + ([] if self.radius_factor < 1.0e-4 else [ lp.GlobalArg("center_int", None, shape=(self.ambient_dim, "nsources"), dim_tags="sep,C"), lp.GlobalArg("center_ext", None, + shape=(self.ambient_dim, "nsources"), dim_tags="sep,C") + ]) + [ + lp.GlobalArg("sources", None, shape=(self.ambient_dim, "nsources"), dim_tags="sep,C"), lp.GlobalArg("proxy_center", None, shape=(self.ambient_dim, "nranges")), @@ -297,8 +307,8 @@ class ProxyGenerator: "..." ], name="find_proxy_radii_knl", - assumptions="dim>=1 and nranges>=1", - fixed_parameters=dict(dim=self.ambient_dim), + assumptions="ndim>=1 and nranges>=1", + fixed_parameters=dict(ndim=self.ambient_dim), lang_version=MOST_RECENT_LANGUAGE_VERSION) knl = lp.tag_inames(knl, "idim*:unr") @@ -313,6 +323,36 @@ class ProxyGenerator: # }}} + def _get_proxy_centers_and_radii(self, actx, source_dd, indices, **kwargs): + from pytential import bind, sym + source_dd = sym.as_dofdesc(source_dd) + discr = self.places.get_discretization( + source_dd.geometry, source_dd.discr_stage) + + from meshmode.dof_array import flatten, thaw + context = { + "sources": flatten(thaw(actx, discr.nodes())), + "srcindices": indices.indices, + "srcranges": indices.ranges + } + + if self.radius_factor > 1.0e-14: + radii = bind(self.places, sym.expansion_radii( + self.ambient_dim, dofdesc=source_dd))(actx) + center_int = bind(self.places, sym.expansion_centers( + self.ambient_dim, -1, dofdesc=source_dd))(actx) + center_ext = bind(self.places, sym.expansion_centers( + self.ambient_dim, +1, dofdesc=source_dd))(actx) + + context["center_int"] = flatten(center_int) + context["center_ext"] = flatten(center_ext) + context["expansion_radii"] = flatten(radii) + + knl = self.get_kernel() + _, (centers, radii,) = knl(actx.queue, **context, **kwargs) + + return centers, radii + def __call__(self, actx, source_dd, indices, **kwargs): """Generate proxy points for each given range of source points in the discretization in *source_dd*. @@ -336,27 +376,8 @@ class ProxyGenerator: # {{{ compute proxy centers and radii - from pytential import bind, sym - source_dd = sym.as_dofdesc(source_dd) - discr = self.places.get_discretization( - source_dd.geometry, source_dd.discr_stage) - - radii = bind(self.places, sym.expansion_radii( - self.ambient_dim, dofdesc=source_dd))(actx) - center_int = bind(self.places, sym.expansion_centers( - self.ambient_dim, -1, dofdesc=source_dd))(actx) - center_ext = bind(self.places, sym.expansion_centers( - self.ambient_dim, +1, dofdesc=source_dd))(actx) - - from meshmode.dof_array import flatten, thaw - knl = self.get_kernel() - _, (centers_dev, radii_dev,) = knl(actx.queue, - sources=flatten(thaw(actx, discr.nodes())), - center_int=flatten(center_int), - center_ext=flatten(center_ext), - expansion_radii=flatten(radii), - srcindices=indices.indices, - srcranges=indices.ranges, **kwargs) + centers_dev, radii_dev = self._get_proxy_centers_and_radii(actx, + source_dd, indices, **kwargs) from pytential.utils import flatten_to_numpy centers = np.vstack(flatten_to_numpy(actx, centers_dev)) @@ -370,7 +391,7 @@ class ProxyGenerator: return np.dot(A, v) + b nproxy = self.nproxy * indices.nblocks - proxies = np.empty((discr.ambient_dim, nproxy), dtype=centers.dtype) + proxies = np.empty((self.ambient_dim, nproxy), dtype=centers.dtype) pxy_nr_base = 0 for i in range(indices.nblocks): @@ -423,13 +444,6 @@ def gather_block_neighbor_points(actx, discr, indices, pxy, indices = indices.get(actx.queue) - # NOTE: this is constructed for multiple reasons: - # * TreeBuilder takes object arrays - # * `srcindices` can be a small subset of nodes, so this will save - # some work - # * `srcindices` may reorder the array returned by nodes(), so this - # makes sure that we have the same order in tree.user_source_ids - # and friends from pytential.utils import flatten_to_numpy sources = flatten_to_numpy(actx, discr.nodes()) sources = make_obj_array([ diff --git a/test/test_linalg_proxy.py b/test/test_linalg_proxy.py index 629e221e..34ed458e 100644 --- a/test/test_linalg_proxy.py +++ b/test/test_linalg_proxy.py @@ -152,7 +152,9 @@ def test_partition_points(ctx_factory, tree_kind, case, visualize=False): @pytest.mark.parametrize("case", PROXY_TEST_CASES) @pytest.mark.parametrize("index_sparsity_factor", [1.0, 0.6]) -def test_proxy_generator(ctx_factory, case, index_sparsity_factor, visualize=False): +@pytest.mark.parametrize("proxy_radius_factor", [0.0, 1.1]) +def test_proxy_generator(ctx_factory, case, + index_sparsity_factor, proxy_radius_factor, visualize=False): """Tests that the proxies generated are all at the correct radius from the points in the cluster, etc. """ @@ -161,7 +163,9 @@ def test_proxy_generator(ctx_factory, case, index_sparsity_factor, visualize=Fal queue = cl.CommandQueue(ctx) actx = PyOpenCLArrayContext(queue) - case = case.copy(index_sparsity_factor=index_sparsity_factor) + case = case.copy( + index_sparsity_factor=index_sparsity_factor, + proxy_radius_factor=proxy_radius_factor) logger.info("\n%s", case) # {{{ check proxies @@ -173,7 +177,10 @@ def test_proxy_generator(ctx_factory, case, index_sparsity_factor, visualize=Fal srcindices = case.get_block_indices(actx, density_discr, matrix_indices=False) from pytential.linalg.proxy import ProxyGenerator - proxies = ProxyGenerator(places)(actx, places.auto_source, srcindices) + proxies = ProxyGenerator(places, + radius_factor=case.proxy_radius_factor, + approx_nproxy=case.proxy_approx_count)( + actx, places.auto_source, srcindices) proxies = proxies.get(queue) pxypoints = np.vstack(proxies.points) @@ -187,10 +194,13 @@ def test_proxy_generator(ctx_factory, case, index_sparsity_factor, visualize=Fal ipxy = proxies.indices.block_indices(i) r = la.norm(pxypoints[:, ipxy] - pxycenters[:, i].reshape(-1, 1), axis=0) - assert np.allclose(r - proxies.radii[i], 0.0, atol=1.0e-14) + p_error = la.norm(r - proxies.radii[i]) r = la.norm(nodes[:, isrc] - pxycenters[:, i].reshape(-1, 1), axis=0) - assert np.all(r - proxies.radii[i] < 0.0) + n_error = la.norm(r - proxies.radii[i], np.inf) + + assert p_error < 1.0e-14, f"block {i}" + assert n_error - proxies.radii[i] < 1.0e-14, f"block {i}" # }}} @@ -216,30 +226,32 @@ def test_proxy_generator(ctx_factory, case, index_sparsity_factor, visualize=Fal bind(places, sym.expansion_radii(ambient_dim))(actx) ) + fig = pt.figure(figsize=(10, 8)) for i in range(srcindices.nblocks): isrc = srcindices.block_indices(i) ipxy = proxies.indices.block_indices(i) - pt.figure(figsize=(10, 8)) - axis = pt.gca() + ax = pt.gca() for j in isrc: c = pt.Circle(ci[:, j], r[j], color="k", alpha=0.1) - axis.add_artist(c) + ax.add_artist(c) c = pt.Circle(ce[:, j], r[j], color="k", alpha=0.1) - axis.add_artist(c) + ax.add_artist(c) - pt.plot(nodes[0], nodes[1], "ko", ms=2.0, alpha=0.5) - pt.plot(nodes[0, srcindices.indices], nodes[1, srcindices.indices], + ax.plot(nodes[0], nodes[1], "ko", ms=2.0, alpha=0.5) + ax.plot(nodes[0, srcindices.indices], nodes[1, srcindices.indices], "o", ms=2.0) - pt.plot(nodes[0, isrc], nodes[1, isrc], "o", ms=2.0) - pt.plot(pxypoints[0, ipxy], pxypoints[1, ipxy], "o", ms=2.0) - pt.axis("equal") - pt.xlim([-1.5, 1.5]) - pt.ylim([-1.5, 1.5]) + ax.plot(nodes[0, isrc], nodes[1, isrc], "o", ms=1.0) + ax.plot(pxypoints[0, ipxy], pxypoints[1, ipxy], "o", ms=1.0) + ax.plot(pxycenters[0, i], pxycenters[1, i], "ko", ms=2.0) + ax.set_aspect("equal") + ax.set_xlim([-1.5, 1.5]) + ax.set_ylim([-1.5, 1.5]) filename = f"test_proxy_generator_{ambient_dim}d_{i:04}" - pt.savefig(filename, dpi=300) - pt.clf() + fig.savefig(filename, dpi=300) + fig.clf() + pt.close(fig) else: from meshmode.discretization.visualization import make_visualizer from meshmode.mesh.processing import ( # noqa @@ -275,7 +287,9 @@ def test_proxy_generator(ctx_factory, case, index_sparsity_factor, visualize=Fal @pytest.mark.parametrize("case", PROXY_TEST_CASES) @pytest.mark.parametrize("index_sparsity_factor", [1.0, 0.6]) -def test_neighbor_points(ctx_factory, case, index_sparsity_factor, visualize=False): +@pytest.mark.parametrize("proxy_radius_factor", [0.0, 1.1]) +def test_neighbor_points(ctx_factory, case, + index_sparsity_factor, proxy_radius_factor, visualize=False): """Test that neighboring points (inside the proxy balls, but outside the current block/cluster) are actually inside. """ @@ -284,7 +298,9 @@ def test_neighbor_points(ctx_factory, case, index_sparsity_factor, visualize=Fal queue = cl.CommandQueue(ctx) actx = PyOpenCLArrayContext(queue) - case = case.copy(index_sparsity_factor=index_sparsity_factor) + case = case.copy( + index_sparsity_factor=index_sparsity_factor, + proxy_radius_factor=proxy_radius_factor) logger.info("\n%s", case) # {{{ check neighboring points @@ -297,7 +313,10 @@ def test_neighbor_points(ctx_factory, case, index_sparsity_factor, visualize=Fal # generate proxy points from pytential.linalg.proxy import ProxyGenerator - proxies = ProxyGenerator(places)(actx, places.auto_source, srcindices) + proxies = ProxyGenerator(places, + radius_factor=case.proxy_radius_factor, + approx_nproxy=case.proxy_approx_count)( + actx, places.auto_source, srcindices) # get neighboring points from pytential.linalg.proxy import gather_block_neighbor_points @@ -334,23 +353,34 @@ def test_neighbor_points(ctx_factory, case, index_sparsity_factor, visualize=Fal except ImportError: return + pxycenters = np.vstack(proxies.centers) + pxyradii = proxies.radii + + fig = pt.figure(figsize=(10, 10)) for i in range(srcindices.nblocks): isrc = srcindices.block_indices(i) inbr = nbrindices.block_indices(i) - pt.figure(figsize=(10, 10)) - pt.plot(nodes[0], nodes[1], "ko", ms=2.0, alpha=0.5) - pt.plot(nodes[0, srcindices.indices], nodes[1, srcindices.indices], + ax = fig.gca() + + ax.plot(nodes[0], nodes[1], "ko", ms=2.0, alpha=0.5) + ax.plot(nodes[0, srcindices.indices], nodes[1, srcindices.indices], "o", ms=2.0) - pt.plot(nodes[0, isrc], nodes[1, isrc], "o", ms=2.0) - pt.plot(nodes[0, inbr], nodes[1, inbr], "o", ms=2.0) - pt.axis("equal") + ax.plot(nodes[0, isrc], nodes[1, isrc], "o", ms=1.0) + ax.plot(nodes[0, inbr], nodes[1, inbr], "o", ms=1.0) + ax.plot(pxycenters[0, i], pxycenters[1, i], "ko", ms=2.0) + ax.axis("equal") + + c = pt.Circle(pxycenters[:, i], pxyradii[i], color="k", alpha=0.1) + ax.add_artist(c) + pt.xlim([-1.5, 1.5]) pt.ylim([-1.5, 1.5]) filename = f"test_area_query_{ambient_dim}d_{i:04}" - pt.savefig(filename, dpi=300) - pt.clf() + fig.savefig(filename, dpi=300) + fig.clf() + pt.close(fig) elif ambient_dim == 3: from meshmode.discretization.visualization import make_visualizer marker = np.empty(density_discr.ndofs) -- GitLab