Skip to content

Add matrix builders to docs #245

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Nov 15, 2024
Merged

Add matrix builders to docs #245

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
1cf7ffb
docs: add matrix builders to docs
alexfikl Nov 15, 2024
35baada
primitives: ignore mypy errors
alexfikl Nov 15, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions doc/symbolic.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,11 @@ Basic objects

.. automodule:: pytential.symbolic.primitives

Matrix Builders
---------------

.. automodule:: pytential.symbolic.matrix

Binding an operator to a discretization
---------------------------------------

Expand Down
155 changes: 139 additions & 16 deletions pytential/symbolic/matrix.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,6 @@
from __future__ import annotations


__copyright__ = """
Copyright (C) 2015 Andreas Kloeckner
Copyright (C) 2018 Alexandru Fikl
Expand All @@ -23,14 +26,45 @@
THE SOFTWARE.
"""

__doc__ = """
.. autoclass:: MatrixBuilderBase
.. autoclass:: ClusterMatrixBuilderBase
:show-inheritance:

Full Matrix Builders
^^^^^^^^^^^^^^^^^^^^

.. autoclass:: MatrixBuilder
:show-inheritance:

.. autoclass:: P2PMatrixBuilder
:show-inheritance:

Cluster Matrix Builders
^^^^^^^^^^^^^^^^^^^^^^^

.. autoclass:: QBXClusterMatrixBuilder
:show-inheritance:

.. autoclass:: P2PClusterMatrixBuilder
:show-inheritance:
"""

from collections.abc import Sequence
from typing import Any

import numpy as np

from sys import intern

from pytools import memoize_method
from arraycontext import flatten, unflatten
from arraycontext import PyOpenCLArrayContext, flatten, unflatten
from meshmode.dof_array import DOFArray
from meshmode.discretization import Discretization

from pytential.collection import GeometryCollection
from pytential.linalg.utils import TargetAndSourceClusterList
from pytential.symbolic.primitives import Variable
from pytential.symbolic.mappers import EvaluationMapperBase


Expand Down Expand Up @@ -70,7 +104,17 @@ def _get_layer_potential_args(actx, places, expr, context=None, include_args=Non
# {{{ base classes for matrix builders

class MatrixBuilderBase(EvaluationMapperBase):
def __init__(self, actx, dep_expr, other_dep_exprs, dep_discr, places, context):
"""
.. automethod:: __init__
"""

def __init__(self,
actx: PyOpenCLArrayContext,
dep_expr: Variable,
other_dep_exprs: Sequence[Variable],
dep_discr: Discretization,
places: GeometryCollection,
context: dict[str, Any]) -> None:
"""
:arg dep_expr: symbolic expression for the input block column
that the builder is evaluating.
Expand Down Expand Up @@ -240,17 +284,25 @@ def map_call(self, expr):

class ClusterMatrixBuilderBase(MatrixBuilderBase):
"""Evaluate individual clusters of a matrix operator, as defined by a
:class:`~pytential.lingla.TargetAndSourceClusterList`.
:class:`~pytential.linalg.TargetAndSourceClusterList`.

Unlike, e.g. :class:`MatrixBuilder`, matrix cluster builders are
significantly reduced in scope. They are basically just meant
to evaluate linear combinations of layer potential operators.
For example, they do not support composition of operators because we
assume that each operator acts directly on the density.

.. automethod:: __init__
"""

def __init__(self, actx, dep_expr, other_dep_exprs,
dep_discr, places, tgt_src_index, context):
def __init__(self,
actx: PyOpenCLArrayContext,
dep_expr: Variable,
other_dep_exprs: Sequence[Variable],
dep_discr: Discretization,
places: GeometryCollection,
tgt_src_index: TargetAndSourceClusterList,
context: dict[str, Any]) -> None:
"""
:arg tgt_src_index: a :class:`~pytential.linalg.TargetAndSourceClusterList`
class describing which clusters are going to be evaluated.
Expand Down Expand Up @@ -313,8 +365,23 @@ class MatrixBuilderDirectResamplerCacheKey:


class MatrixBuilder(MatrixBuilderBase):
def __init__(self, actx, dep_expr, other_dep_exprs,
dep_discr, places, context, _weighted=True):
"""A matrix builder that evaluates the full QBX-mediated operator.

.. automethod:: __init__
"""

def __init__(self,
actx: PyOpenCLArrayContext,
dep_expr: Variable,
other_dep_exprs: Sequence[Variable],
dep_discr: Discretization,
places: GeometryCollection,
context: dict[str, Any],
_weighted: bool = True) -> None:
"""
:arg _weighted: if *True*, the quadrature weights and area elements are
added to the operator matrix evaluation.
"""
super().__init__(actx, dep_expr, other_dep_exprs, dep_discr, places, context)
self.weighted = _weighted

Expand Down Expand Up @@ -445,9 +512,28 @@ def map_int_g(self, expr):
# {{{ p2p matrix builder

class P2PMatrixBuilder(MatrixBuilderBase):
def __init__(self, actx, dep_expr, other_dep_exprs,
dep_discr, places, context,
exclude_self=True, _weighted=False):
"""A matrix builder that evaluates the full point-to-point kernel interactions.

.. automethod:: __init__
"""

def __init__(self,
actx: PyOpenCLArrayContext,
dep_expr: Variable,
other_dep_exprs: Sequence[Variable],
dep_discr: Discretization,
places: GeometryCollection,
context: dict[str, Any],
exclude_self: bool = True,
_weighted: bool = False) -> None:
"""
:arg exclude_self: if *True*, interactions where the source and target
indices are the same are ignored. This should be *True* in most
cases where this is possible, since the kernels are singular at those
points.
:arg _weighted: if *True*, the quadrature weights and area elements are
added to the operator matrix evaluation.
"""
super().__init__(actx, dep_expr, other_dep_exprs, dep_discr, places, context)

self.exclude_self = exclude_self
Expand Down Expand Up @@ -529,9 +615,26 @@ def map_int_g(self, expr):
# {{{ cluster matrix builders

class QBXClusterMatrixBuilder(ClusterMatrixBuilderBase):
def __init__(self, actx, dep_expr, other_dep_exprs, dep_discr,
places, tgt_src_index, context,
exclude_self=False, _weighted=True):
"""A matrix builder that evaluates a cluster (block) using the QBX method.

.. automethod:: __init__
"""

def __init__(self,
actx: PyOpenCLArrayContext,
dep_expr: Variable,
other_dep_exprs: Sequence[Variable],
dep_discr: Discretization,
places: GeometryCollection,
tgt_src_index: TargetAndSourceClusterList,
context: dict[str, Any],
exclude_self: bool = False,
_weighted: bool = True) -> None:
"""
:arg exclude_self: this argument is ignored.
:arg _weighted: if *True*, the quadrature weights and area elements are
added to the operator matrix evaluation.
"""
super().__init__(
actx, dep_expr, other_dep_exprs, dep_discr, places,
tgt_src_index, context)
Expand Down Expand Up @@ -620,9 +723,29 @@ def map_int_g(self, expr):


class P2PClusterMatrixBuilder(ClusterMatrixBuilderBase):
def __init__(self, actx, dep_expr, other_dep_exprs, dep_discr,
places, tgt_src_index, context,
exclude_self=False, _weighted=False):
"""A matrix builder that evaluates a cluster (block) point-to-point.

.. automethod:: __init__
"""

def __init__(self,
actx: PyOpenCLArrayContext,
dep_expr: Variable,
other_dep_exprs: Sequence[Variable],
dep_discr: Discretization,
places: GeometryCollection,
tgt_src_index: TargetAndSourceClusterList,
context: dict[str, Any],
exclude_self: bool = False,
_weighted: bool = False) -> None:
"""
:arg exclude_self: if *True*, interactions where the source and target
indices are the same are ignored. This should be *True* in most
cases where this is possible, since the kernels are singular at those
points.
:arg _weighted: if *True*, the quadrature weights and area elements are
added to the operator matrix evaluation.
"""
super().__init__(
actx, dep_expr, other_dep_exprs, dep_discr, places,
tgt_src_index, context)
Expand Down
13 changes: 8 additions & 5 deletions pytential/symbolic/primitives.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -558,7 +558,8 @@ class NumReferenceDerivative(DiscretizationProperty):
def __new__(cls,
ref_axes: int | tuple[tuple[int, int], ...] | None = None,
operand: Operand | None = None,
dofdesc: DOFDescriptor | None = None) -> "NumReferenceDerivative":
dofdesc: DOFDescriptor | None = None,
) -> "NumReferenceDerivative":
# If the constructor is handed a multivector object, return an
# object array of the operator applied to each of the
# coefficients in the multivector.
Expand All @@ -567,7 +568,9 @@ def __new__(cls,
def make_op(operand_i):
return cls(ref_axes, operand_i, as_dofdesc(dofdesc))

return componentwise(make_op, operand)
# FIXME: mypy is right: new should return `cls` instances and we're
# abusing it to vectorize the call like this.
return componentwise(make_op, operand) # type: ignore[return-value]
Comment on lines -570 to +573
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not quite sure what else to do about this?

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

https://docs.python.org/3/reference/datamodel.html#object.__new__ says:

If __new__() does not return an instance of cls, then the new instance’s __init__() method will not be invoked.

I read this as implicit permission for __new__ to return whatever it wants.

What happens if you're honest in the type annotation for __new__?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

At least mypy is unhappy

pytential/symbolic/primitives.py:558: error: "__new__" must return a class instance (got "NumReferenceDerivative | ndarray[Any, dtype[Any]]")  [misc]

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also tried doing an overload like

    @overload
    def __new__(cls,
                ref_axes: int | tuple[tuple[int, int], ...] | None = None,
                operand: Expression | None = None,
                dofdesc: DOFDescriptor | None = None,
                ) -> "NumReferenceDerivative": ...

    @overload
    def __new__(cls,
                ref_axes: int | tuple[tuple[int, int], ...] | None = None,
                operand: Union["np.ndarray[Any, np.dtype[Any]]", None] = None,
                dofdesc: DOFDescriptor | None = None,
                ) -> "np.ndarray[Any, np.dtype[Any]]": ...

and it still gives

pytential/symbolic/primitives.py:566: error: Incompatible return type for "__new__" (returns "ndarray[Any, dtype[Any]]", but must return a subtype of "NumReferenceDerivative")  [misc]

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Mypy issue/PR rabbit hole in reverse chronological order:

Including a BDFL quote:

Honestly I don't think that A() should ever return something that's not an instance of A (or a subclass), even though technically new() can return anything it likes.

I don't think there's any way that mypy will learn what's going on here in the short term. Stylistically, it may be better to move away from overloading __new__ and maybe instead have the user call componentwise directly. But that's a different discussion, #246.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

He does agree that this sort of symbolic simplification usage is quite pleasant though 😁
python/mypy#1020 (comment)

else:
return DiscretizationProperty.__new__(cls)

Expand Down Expand Up @@ -1140,7 +1143,7 @@ def __new__(cls,
def make_op(operand_i):
return cls(from_dd, to_dd, operand_i)

return componentwise(make_op, operand)
return componentwise(make_op, operand) # type: ignore[return-value]
else:
return Expression.__new__(cls)

Expand Down Expand Up @@ -1185,7 +1188,7 @@ def __new__(cls,
def make_op(operand_i):
return cls(operand_i)

return componentwise(make_op, operand)
return componentwise(make_op, operand) # type: ignore[return-value]
else:
return Expression.__new__(cls)

Expand Down Expand Up @@ -1248,7 +1251,7 @@ def __new__(cls,
def make_op(operand_i):
return cls(operand_i, as_dofdesc(dofdesc))

return componentwise(make_op, operand)
return componentwise(make_op, operand) # type: ignore[return-value]
else:
return Expression.__new__(cls)

Expand Down
Loading