From 38ecebd44b1891b6143eaa9c04525e78e3574322 Mon Sep 17 00:00:00 2001 From: Ben Volokh <89551265+benv123@users.noreply.github.com> Date: Fri, 14 Jul 2023 07:22:44 -0400 Subject: [PATCH] Simplex Domain (#152) * add Simplex Domain * unittests for Simplex * docs update --------- Co-authored-by: Dario Coscia <93731561+dario-coscia@users.noreply.github.com> --- pina/geometry/__init__.py | 4 +- pina/geometry/simplex.py | 225 ++++++++++++++++++++++++++++ tests/test_geometry/test_simplex.py | 141 +++++++++++++++++ 3 files changed, 369 insertions(+), 1 deletion(-) create mode 100644 pina/geometry/simplex.py create mode 100644 tests/test_geometry/test_simplex.py diff --git a/pina/geometry/__init__.py b/pina/geometry/__init__.py index 61253db..5d38b32 100644 --- a/pina/geometry/__init__.py +++ b/pina/geometry/__init__.py @@ -6,7 +6,8 @@ __all__ = [ 'Intersection', 'Exclusion', 'Difference', - 'OperationInterface' + 'OperationInterface', + 'SimplexDomain' ] from .location import Location @@ -17,3 +18,4 @@ from .intersection_domain import Intersection from .union_domain import Union from .difference_domain import Difference from .operation_interface import OperationInterface +from .simplex import SimplexDomain diff --git a/pina/geometry/simplex.py b/pina/geometry/simplex.py new file mode 100644 index 0000000..e3670a2 --- /dev/null +++ b/pina/geometry/simplex.py @@ -0,0 +1,225 @@ +import torch +from .location import Location +from pina.geometry import CartesianDomain +from pina import LabelTensor +from ..utils import check_consistency + + +class SimplexDomain(Location): + """PINA implementation of a Simplex.""" + + def __init__(self, simplex_matrix, sample_surface=False): + """ + :param simplex_matrix: A matrix of LabelTensor objects representing + a vertex of the simplex (a tensor), and the coordinates of the + point (a list of labels). + :type simplex_matrix: list[LabelTensor] + :param sample_surface: A variable for choosing sample strategies. If + `sample_surface=True` only samples on the Simplex surface + frontier are taken. If `sample_surface=False`, no such criteria + is followed. + :type sample_surface: bool + .. warning:: + Sampling for dimensions greater or equal to 10 could result + in a shrinking of the simplex, which degrades the quality + of the samples. For dimensions higher than 10, other algorithms + for sampling should be used. + :Example: + >>> spatial_domain = SimplexDomain( + [ + LabelTensor(torch.tensor([[0, 0]]), labels=["x", "y"]), + LabelTensor(torch.tensor([[1, 1]]), labels=["x", "y"]), + LabelTensor(torch.tensor([[0, 2]]), labels=["x", "y"]), + ], sample_surface = True + ) + """ + + # check consistency of sample_surface as bool + check_consistency(sample_surface, bool) + self._sample_surface = sample_surface + + # check consistency of simplex_matrix as list or tuple + check_consistency([simplex_matrix], (list, tuple)) + + # check everything within simplex_matrix is a LabelTensor + check_consistency(simplex_matrix, LabelTensor) + + # check consistency of labels + self._coordinates = simplex_matrix[0].labels + if not all(vertex.labels == self._coordinates for vertex in simplex_matrix): + raise ValueError(f"Labels don't match.") + + # creating vertices matrix + self._vertices_matrix = torch.cat(simplex_matrix) + self._vertices_matrix.labels = self._coordinates + + # creating basis vectors for simplex + self._vectors_shifted = ( + (self._vertices_matrix.T)[:, :-1] - (self._vertices_matrix.T)[:, None, -1] + ) + + # build cartesian_bound + self._cartesian_bound = self._build_cartesian(self.variables) + + @property + def coordinates(self): + return self._coordinates + + @property + def variables(self): + return self._vertices_matrix + + def _build_cartesian(self, vertices): + """ + Build Cartesian border for Simplex domain to be used in sampling. + :param vertex_matrix: matrix of vertices + :type vertices: list[list] + :return: Cartesian border for triangular domain + :rtype: CartesianDomain + """ + + span_dict = {} + + for i, coord in enumerate(self._coordinates): + sorted_vertices = sorted(vertices, key=lambda vertex: vertex[i]) + # respective coord bounded by the lowest and highest values + span_dict[coord] = [sorted_vertices[0][i], sorted_vertices[-1][i]] + + return CartesianDomain(span_dict) + + def is_inside(self, point, check_border=False): + """ + Check if a point is inside the simplex. + Uses the algorithm described involving barycentric coordinates: + https://en.wikipedia.org/wiki/Barycentric_coordinate_system. + .. note:: + When ```'sample_surface'``` in the ```'__init()__'``` + is set to ```'True'```, then the method only checks + points on the surface, and not inside the domain. + :param point: Point to be checked. + :type point: LabelTensor + :param check_border: Check if the point is also on the frontier + of the simplex, default False. + :type check_border: bool + :return: Returning True if the point is inside, False otherwise. + :rtype: bool + """ + + if not all([label in self.coordinates for label in point.labels]): + raise ValueError( + "Point labels different from constructor" + f" dictionary labels. Got {point.labels}," + f" expected {self.variables}." + ) + + # shift point + point_shift = point.T - (self._vertices_matrix.T)[:, None, -1] + + # compute barycentric coordinates + lambda_ = torch.linalg.solve(self._vectors_shifted * 1.0, point_shift * 1.0) + lambda_1 = 1.0 - torch.sum(lambda_) + lambdas = torch.vstack([lambda_, lambda_1]) + + # perform checks + if not check_border: + return all(torch.gt(lambdas, 0.0)) and all(torch.lt(lambdas, 1.0)) + + return all(torch.ge(lambdas, 0)) and ( + any(torch.eq(lambdas, 0)) or any(torch.eq(lambdas, 1)) + ) + + def _sample_interior_randomly(self, n, variables): + """ + Randomly sample points inside a simplex of arbitrary + dimension, without the boundary. + :param int n: Number of points to sample in the shape. + :param variables: pinn variable to be sampled, defaults to 'all'. + :type variables: str or list[str], optional + :return: Returns tensor of n sampled points. + :rtype: torch.Tensor + """ + + # =============== For Developers ================ # + # + # The sampling startegy used is fairly simple. + # First we sample a random vector from the hypercube + # which contains the simplex. Then, if the point + # sampled is inside the simplex, we add it as a valid + # one. + # + # =============================================== # + + sampled_points = [] + while len(sampled_points) < n: + sampled_point = self._cartesian_bound.sample( + n=1, mode="random", variables=variables + ) + + if self.is_inside(sampled_point, self._sample_surface): + sampled_points.append(sampled_point) + return torch.cat(sampled_points, dim=0) + + def _sample_boundary_randomly(self, n): + """ + Randomly sample points on the boundary of a simplex + of arbitrary dimensions. + :param int n: Number of points to sample in the shape. + :return: Returns tensor of n sampled points + :rtype: torch.Tensor + """ + + # =============== For Developers ================ # + # + # The sampling startegy used is fairly simple. + # We first sample the lambdas in [0, 1] domain, + # we then set to zero only one lambda, and normalize. + # Finally, we compute the matrix product between the + # lamdas and the vertices matrix. + # + # =============================================== # + + sampled_points = [] + + while len(sampled_points) < n: + # extract number of vertices + number_of_vertices = self._vertices_matrix.shape[1] + # extract idx lambda to set to zero randomly + idx_lambda = torch.randint(low=0, high=number_of_vertices, size=(1,)) + # build lambda vector + # 1. sampling [1, 2) + lambdas = torch.rand((number_of_vertices, 1)) + # 2. setting lambdas[idx_lambda] to 0 + lambdas[idx_lambda] = 0 + # 3. normalize + lambdas /= lambdas.sum() + # 4. compute dot product + sampled_points.append(self._vertices_matrix @ lambdas) + + return torch.cat(sampled_points, dim=1).T + + def sample(self, n, mode="random", variables="all"): + """ + Sample n points from Simplex domain. + :param int n: Number of points to sample in the shape. + :param str mode: Mode for sampling, defaults to 'random'. + Available modes include: 'random'. + :param variables: pinn variable to be sampled, defaults to 'all'. + :type variables: str or list[str], optional + :return: Returns LabelTensor of n sampled points + :rtype: LabelTensor(tensor) + .. warning:: + When ``sample_surface = True`` in the initialization, all + the variables are sampled, despite passing different once + in ``variables``. + """ + + if mode in ["random"]: + if self._sample_surface: + sample_pts = self._sample_boundary_randomly(n) + else: + sample_pts = self._sample_interior_randomly(n, variables) + + else: + raise NotImplementedError(f"mode={mode} is not implemented.") + + return LabelTensor(sample_pts, labels=self.variables) diff --git a/tests/test_geometry/test_simplex.py b/tests/test_geometry/test_simplex.py new file mode 100644 index 0000000..e0cf081 --- /dev/null +++ b/tests/test_geometry/test_simplex.py @@ -0,0 +1,141 @@ +import torch +import pytest + +from pina import LabelTensor +from pina.geometry import SimplexDomain + + +def test_constructor(): + SimplexDomain( + [ + LabelTensor(torch.tensor([[0, 0]]), labels=["x", "y"]), + LabelTensor(torch.tensor([[1, 1]]), labels=["x", "y"]), + LabelTensor(torch.tensor([[0, 2]]), labels=["x", "y"]), + ] + ) + SimplexDomain( + [ + LabelTensor(torch.tensor([[0, 0]]), labels=["x", "y"]), + LabelTensor(torch.tensor([[1, 1]]), labels=["x", "y"]), + LabelTensor(torch.tensor([[0, 2]]), labels=["x", "y"]), + ], + sample_surface=True, + ) + with pytest.raises(ValueError): + SimplexDomain( + [ + LabelTensor(torch.tensor([[0, 0]]), labels=["x", "y"]), + LabelTensor(torch.tensor([[1, 1]]), labels=["x", "z"]), + LabelTensor(torch.tensor([[0, 2]]), labels=["x", "a"]), + ] + ) + SimplexDomain( + [ + LabelTensor(torch.tensor([[0, 0]]), labels=["x", "y"]), + [1, 1], + LabelTensor(torch.tensor([[0, 2]]), labels=["x", "y"]), + ] + ) + + +def test_is_inside_faulty_point(): + domain = SimplexDomain( + [ + LabelTensor(torch.tensor([[0, 0]]), labels=["x", "y"]), + LabelTensor(torch.tensor([[2, 2]]), labels=["x", "y"]), + LabelTensor(torch.tensor([[4, 0]]), labels=["x", "y"]), + ], + sample_surface=True, + ) + pt = LabelTensor(torch.tensor([[2, 1]]), ["x", "z"]) + with pytest.raises(ValueError): + assert domain.is_inside(point=pt, check_border=False) == True + + +def test_is_inside_2D_check_border_true(): + domain = SimplexDomain( + [ + LabelTensor(torch.tensor([[0, 0]]), labels=["x", "y"]), + LabelTensor(torch.tensor([[2, 2]]), labels=["x", "y"]), + LabelTensor(torch.tensor([[4, 0]]), labels=["x", "y"]), + ], + sample_surface=True, + ) + pt1 = LabelTensor(torch.tensor([[0, 0]]), ["x", "y"]) + pt2 = LabelTensor(torch.tensor([[2, 2]]), ["x", "y"]) + pt3 = LabelTensor(torch.tensor([[4, 0]]), ["x", "y"]) + pt4 = LabelTensor(torch.tensor([[3, 1]]), ["x", "y"]) + pt5 = LabelTensor(torch.tensor([[2, 1]]), ["x", "y"]) + pt6 = LabelTensor(torch.tensor([[100, 100]]), ["x", "y"]) + pts = [pt1, pt2, pt3, pt4, pt5, pt6] + for pt, exp_result in zip(pts, [True, True, True, True, False, False]): + assert domain.is_inside(point=pt, check_border=True) == exp_result + + +def test_is_inside_2D_check_border_false(): + domain = SimplexDomain( + [ + LabelTensor(torch.tensor([[0, 0]]), labels=["x", "y"]), + LabelTensor(torch.tensor([[2, 2]]), labels=["x", "y"]), + LabelTensor(torch.tensor([[4, 0]]), labels=["x", "y"]), + ], + sample_surface=False, + ) + pt1 = LabelTensor(torch.tensor([[0, 0]]), ["x", "y"]) + pt2 = LabelTensor(torch.tensor([[2, 2]]), ["x", "y"]) + pt3 = LabelTensor(torch.tensor([[4, 0]]), ["x", "y"]) + pt4 = LabelTensor(torch.tensor([[3, 1]]), ["x", "y"]) + pt5 = LabelTensor(torch.tensor([[2, 1]]), ["x", "y"]) + pt6 = LabelTensor(torch.tensor([[2.5, 1]]), ["x", "y"]) + pt7 = LabelTensor(torch.tensor([[100, 100]]), ["x", "y"]) + pts = [pt1, pt2, pt3, pt4, pt5, pt6, pt7] + for pt, exp_result in zip(pts, [False, False, False, False, True, True, False]): + assert domain.is_inside(point=pt, check_border=False) == exp_result + + +def test_is_inside_3D_check_border_true(): + domain = SimplexDomain( + [ + LabelTensor(torch.tensor([[0, 0, 0]]), labels=["x", "y", "z"]), + LabelTensor(torch.tensor([[2, 2, 0]]), labels=["x", "y", "z"]), + LabelTensor(torch.tensor([[4, 0, 0]]), labels=["x", "y", "z"]), + LabelTensor(torch.tensor([[0, 0, 20]]), labels=["x", "y", "z"]), + ], + sample_surface=True, + ) + pt1 = LabelTensor(torch.tensor([[0, 0, 0]]), ["x", "y", "z"]) + pt2 = LabelTensor(torch.tensor([[2, 2, 0]]), ["x", "y", "z"]) + pt3 = LabelTensor(torch.tensor([[4, 0, 0]]), ["x", "y", "z"]) + pt4 = LabelTensor(torch.tensor([[3, 1, 0]]), ["x", "y", "z"]) + pt5 = LabelTensor(torch.tensor([[2, 1, 0]]), ["x", "y", "z"]) + pt6 = LabelTensor(torch.tensor([[100, 100, 0]]), ["x", "y", "z"]) + pt7 = LabelTensor(torch.tensor([[0, 0, 19]]), ["x", "y", "z"]) + pt8 = LabelTensor(torch.tensor([[0, 0, 20]]), ["x", "y", "z"]) + pt9 = LabelTensor(torch.tensor([[2, 1, 1]]), ["x", "y", "z"]) + pts = [pt1, pt2, pt3, pt4, pt5, pt6, pt7, pt8, pt9] + for pt, exp_result in zip( + pts, [True, True, True, True, True, False, True, True, False] + ): + assert domain.is_inside(point=pt, check_border=True) == exp_result + + +def test_is_inside_3D_check_border_false(): + domain = SimplexDomain( + [ + LabelTensor(torch.tensor([[0, 0, 0]]), labels=["x", "y", "z"]), + LabelTensor(torch.tensor([[2, 2, 0]]), labels=["x", "y", "z"]), + LabelTensor(torch.tensor([[4, 0, 0]]), labels=["x", "y", "z"]), + LabelTensor(torch.tensor([[0, 0, 20]]), labels=["x", "y", "z"]), + ], + sample_surface=False, + ) + pt1 = LabelTensor(torch.tensor([[0, 0, 0]]), ["x", "y", "z"]) + pt2 = LabelTensor(torch.tensor([[3, 1, 0]]), ["x", "y", "z"]) + pt3 = LabelTensor(torch.tensor([[2, 1, 0]]), ["x", "y", "z"]) + pt4 = LabelTensor(torch.tensor([[100, 100, 0]]), ["x", "y", "z"]) + pt5 = LabelTensor(torch.tensor([[0, 0, 19]]), ["x", "y", "z"]) + pt6 = LabelTensor(torch.tensor([[0, 0, 20]]), ["x", "y", "z"]) + pt7 = LabelTensor(torch.tensor([[2, 1, 1]]), ["x", "y", "z"]) + pts = [pt1, pt2, pt3, pt4, pt5, pt6, pt7] + for pt, exp_result in zip(pts, [False, False, False, False, False, False, True]): + assert domain.is_inside(point=pt, check_border=False) == exp_result