Solving problems related to Geometry (#118)

* fix and add tests
* minor fix on domain classes

---------

Co-authored-by: Dario Coscia <dariocoscia@cli-10-110-0-208.WIFIeduroamSTUD.units.it>
Co-authored-by: Dario Coscia <dariocoscia@dhcp-040.eduroam.sissa.it>

pina/geometry/cartesian.py

@@ -247,23 +247,27 @@ class CartesianDomain(Location):
         :param point: Point to be checked
         :type point: LabelTensor
         :param check_border: Check if the point is also on the frontier
-            of the ellipsoid, default False.
+            of the hypercube, default False.
         :type check_border: bool
         :return: Returning True if the point is inside, False otherwise.
         :rtype: bool
         """
         is_inside = []
+
+        # check fixed variables
+        for variable, value in self.fixed_.items():
+            if variable in point.labels:
+                is_inside.append(point.extract([variable]) == value)
+
+        # check not fixed variables
         for variable, bound in self.range_.items():
             if variable in point.labels:
-                if bound[0] <= point.extract([variable]) <= bound[1]:
-                    is_inside.append(True)
+
+                if check_border:
+                    check = bound[0] <= point.extract([variable]) <= bound[1]
                 else:
-                    is_inside.append(False)
+                    check = bound[0] < point.extract([variable]) < bound[1]
+            
+                is_inside.append(check)
 
         return all(is_inside)
-
-        # TODO check the fixed_ dimensions
-        # for variable, value in self.fixed_.items():
-        #     if variable in point.labels:
-        #         if not (point.extract[variable] == value):
-        #             return False

pina/geometry/difference_domain.py

@@ -16,13 +16,13 @@ class Difference(Location):
     def sample(self, n, mode='random', variables='all'):
         """
         """
-        assert mode is 'random', 'Only random mode is implemented'
+        assert mode == 'random', 'Only random mode is implemented'
 
         samples = []
         while len(samples) < n:
             sample = self.first.sample(1, 'random')
             if not self.second.is_inside(sample):
-                samples.append(sample.tolist()[0])
+                samples.append(sample)
 
         import torch
-        return LabelTensor(torch.tensor(samples), labels=['x', 'y'])
+        return LabelTensor(torch.cat(samples), labels=['x', 'y'])

pina/geometry/ellipsoid.py

@@ -2,6 +2,7 @@ import torch
 
 from .location import Location
 from ..label_tensor import LabelTensor
+from ..utils import check_consistency
 
 
 class EllipsoidDomain(Location):
@@ -39,9 +40,8 @@ class EllipsoidDomain(Location):
         self._centers = None
         self._axis = None
 
-        if not isinstance(sample_surface, bool):
-            raise ValueError('sample_surface must be bool type.')
-
+        # checking consistency
+        check_consistency(sample_surface, bool)
         self._sample_surface = sample_surface
 
         for k, v in ellipsoid_dict.items():
@@ -81,9 +81,14 @@ class EllipsoidDomain(Location):
         return list(self.fixed_.keys()) + list(self.range_.keys())
 
     def is_inside(self, point, check_border=False):
-        """Check if a point is inside the ellipsoid.
+        """Check if a point is inside the ellipsoid domain.
 
-        :param point: Point to be checked
+        .. 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 ellipsoid, default False.
@@ -92,29 +97,40 @@ class EllipsoidDomain(Location):
         :rtype: bool
         """
 
-        if not isinstance(point, LabelTensor):
-            raise ValueError('point expected to be LabelTensor.')
-
-        # get axis ellipse
+        # small check that point is labeltensor
+        check_consistency(point, LabelTensor)
+        
+        # get axis ellipse as tensors
         list_dict_vals = list(self._axis.values())
         tmp = torch.tensor(list_dict_vals, dtype=torch.float)
-        ax_sq = LabelTensor(tmp.reshape(1, -1)**2, list(self._axis.keys()))
+        ax_sq = LabelTensor(tmp.reshape(1, -1)**2, self.variables)
+
+        # get centers ellipse as tensors
+        list_dict_vals = list(self._centers.values())
+        tmp = torch.tensor(list_dict_vals, dtype=torch.float)
+        centers = LabelTensor(tmp.reshape(1, -1), self.variables)
 
         if not all([i in ax_sq.labels for i in point.labels]):
             raise ValueError('point labels different from constructor'
                              f' dictionary labels. Got {point.labels},'
                              f' expected {ax_sq.labels}.')
 
-        # point square
-        point_sq = point.pow(2)
+        # point square + shift center
+        point_sq = (point - centers).pow(2)
         point_sq.labels = point.labels
 
         # calculate ellispoid equation
         eqn = torch.sum(point_sq.extract(ax_sq.labels) / ax_sq) - 1.
 
+        # if we have sampled only the surface, we check that the 
+        # point is inside the surface border only
+        if self._sample_surface:
+            return torch.allclose(eqn, torch.zeros_like(eqn))
+
+        # otherwise we check the ellipse 
         if check_border:
             return bool(eqn <= 0)
-
+        
         return bool(eqn < 0)
 
     def _sample_range(self, n, mode, variables):
@@ -265,4 +281,4 @@ class EllipsoidDomain(Location):
         if mode in ['random']:
             return _Nd_sampler(n, mode, variables)
         else:
-            raise ValueError(f'mode={mode} is not valid.')
+            raise NotImplemented(f'mode={mode} is not implemented.')

pina/geometry/union_domain.py

@@ -2,6 +2,7 @@ import torch
 from .location import Location
 from ..utils import check_consistency
 from ..label_tensor import LabelTensor
+import random
 
 
 class Union(Location):
@@ -87,7 +88,7 @@ class Union(Location):
             >>> ellipsoid2 = EllipsoidDomain({'x': [0, 2], 'y': [0, 2]})
 
             # Create a union of the ellipsoid domains
-            >>> union = GeometryUnion([ellipsoid1, ellipsoid2])
+            >>> union = Union([ellipsoid1, ellipsoid2])
 
             >>> union.sample(n=1000)
                 LabelTensor([[-0.2025,  0.0072],
@@ -108,11 +109,18 @@ class Union(Location):
         num_points = n // len(self.geometries)
 
         # sample the points
-        for i, geometry in enumerate(self.geometries):
-            # add to sample total if remainder is not 0
-            if i < remainder:
-                num_points += 1
-            sampled_points.append(geometry.sample(num_points, mode, variables))
+        # NB. geometries as shuffled since if we sample
+        # multiple times just one point, we would end
+        # up sampling only from the first geometry.
+        iter_ = random.sample(self.geometries, len(self.geometries))
+        for i, geometry in enumerate(iter_):
+            # int(i < remainder) is one only if we have a remainder
+            # different than zero. Notice that len(geometries) is
+            # always smaller than remaider.
+            sampled_points.append(geometry.sample(num_points + int(i < remainder), mode, variables))
+            # in case number of sampled points is smaller than the number of geometries
+            if len(sampled_points) >= n:
+                break
 
         return LabelTensor(torch.cat(sampled_points), labels=[f'{i}' for i in self.variables])