equation class, fix minor bugs, diff domain (#89)

* equation class
* difference domain
* dummy dataloader
* writer class
* refactoring and minor fix
This commit is contained in:
Nicola Demo
2023-05-15 16:06:01 +02:00
parent be11110bb2
commit 0e3625de80
25 changed files with 691 additions and 246 deletions

10
pina/geometry/__init__.py Normal file
View File

@@ -0,0 +1,10 @@
__all__ = [
'Location',
'CartesianDomain',
'EllipsoidDomain',
]
from .location import Location
from .cartesian import CartesianDomain
from .ellipsoid import EllipsoidDomain
from .difference_domain import Difference

269
pina/geometry/cartesian.py Normal file
View File

@@ -0,0 +1,269 @@
import torch
from .location import Location
from ..label_tensor import LabelTensor
from ..utils import torch_lhs, chebyshev_roots
class CartesianDomain(Location):
"""PINA implementation of Hypercube domain."""
def __init__(self, span_dict):
"""
:param span_dict: A dictionary with dict-key a string representing
the input variables for the pinn, and dict-value a list with
the domain extrema.
:type span_dict: dict
:Example:
>>> spatial_domain = Span({'x': [0, 1], 'y': [0, 1]})
"""
self.fixed_ = {}
self.range_ = {}
for k, v in span_dict.items():
if isinstance(v, (int, float)):
self.fixed_[k] = v
elif isinstance(v, (list, tuple)) and len(v) == 2:
self.range_[k] = v
else:
raise TypeError
@property
def variables(self):
"""Spatial variables.
:return: Spatial variables defined in '__init__()'
:rtype: list[str]
"""
return list(self.fixed_.keys()) + list(self.range_.keys())
def update(self, new_span):
"""Adding new dimensions on the span
:param new_span: A new span object to merge
:type new_span: Span
:Example:
>>> spatial_domain = Span({'x': [0, 1], 'y': [0, 1]})
>>> spatial_domain.variables
['x', 'y']
>>> spatial_domain_2 = Span({'z': [3, 4], 'w': [0, 1]})
>>> spatial_domain.update(spatial_domain_2)
>>> spatial_domain.variables
['x', 'y', 'z', 'w']
"""
self.fixed_.update(new_span.fixed_)
self.range_.update(new_span.range_)
def _sample_range(self, n, mode, bounds):
"""Rescale the samples to the correct bounds
:param n: Number of points to sample, see Note below
for reference.
:type n: int
:param mode: Mode for sampling, defaults to 'random'.
Available modes include: random sampling, 'random';
latin hypercube sampling, 'latin' or 'lh';
chebyshev sampling, 'chebyshev'; grid sampling 'grid'.
:type mode: str, optional
:param bounds: Bounds to rescale the samples.
:type bounds: torch.Tensor
:return: Rescaled sample points.
:rtype: torch.Tensor
"""
dim = bounds.shape[0]
if mode in ['chebyshev', 'grid'] and dim != 1:
raise RuntimeError('Something wrong in Span...')
if mode == 'random':
pts = torch.rand(size=(n, dim))
elif mode == 'chebyshev':
pts = chebyshev_roots(n).mul(.5).add(.5).reshape(-1, 1)
elif mode == 'grid':
pts = torch.linspace(0, 1, n).reshape(-1, 1)
elif mode == 'lh' or mode == 'latin':
pts = torch_lhs(n, dim)
pts *= bounds[:, 1] - bounds[:, 0]
pts += bounds[:, 0]
return pts
def sample(self, n, mode='random', variables='all'):
"""Sample routine.
:param n: Number of points to sample, see Note below
for reference.
:type n: int
:param mode: Mode for sampling, defaults to 'random'.
Available modes include: random sampling, 'random';
latin hypercube sampling, 'latin' or 'lh';
chebyshev sampling, 'chebyshev'; grid sampling 'grid'.
:type mode: str, optional
:param variables: pinn variable to be sampled, defaults to 'all'.
:type variables: str or list[str], optional
.. note::
The total number of points sampled in case of multiple variables
is not 'n', and it depends on the chosen 'mode'. If 'mode' is
'grid' or 'chebyshev', the points are sampled independentely
across the variables and the results crossed together, i.e. the
final number of points is 'n' to the power of the number of
variables. If 'mode' is 'random', 'lh' or 'latin', the variables
are sampled all together, and the final number of points
.. warning::
The extrema values of Span are always sampled only for 'grid' mode.
:Example:
>>> spatial_domain = Span({'x': [0, 1], 'y': [0, 1]})
>>> spatial_domain.sample(n=4, mode='random')
tensor([[0.0108, 0.7643],
[0.4477, 0.8015],
[0.2063, 0.8087],
[0.8735, 0.6349]])
>>> spatial_domain.sample(n=4, mode='grid')
tensor([[0.0000, 0.0000],
[0.3333, 0.0000],
[0.6667, 0.0000],
[1.0000, 0.0000],
[0.0000, 0.3333],
[0.3333, 0.3333],
[0.6667, 0.3333],
[1.0000, 0.3333],
[0.0000, 0.6667],
[0.3333, 0.6667],
[0.6667, 0.6667],
[1.0000, 0.6667],
[0.0000, 1.0000],
[0.3333, 1.0000],
[0.6667, 1.0000],
[1.0000, 1.0000]])
"""
def _1d_sampler(n, mode, variables):
""" Sample independentely the variables and cross the results"""
tmp = []
for variable in variables:
if variable in self.range_.keys():
bound = torch.tensor([self.range_[variable]])
pts_variable = self._sample_range(n, mode, bound)
pts_variable = pts_variable.as_subclass(LabelTensor)
pts_variable.labels = [variable]
tmp.append(pts_variable)
result = tmp[0]
for i in tmp[1:]:
result = result.append(i, mode='cross')
for variable in variables:
if variable in self.fixed_.keys():
value = self.fixed_[variable]
pts_variable = torch.tensor([[value]]).repeat(
result.shape[0], 1)
pts_variable = pts_variable.as_subclass(LabelTensor)
pts_variable.labels = [variable]
result = result.append(pts_variable, mode='std')
return result
def _Nd_sampler(n, mode, variables):
"""Sample all the variables together
:param n: Number of points to sample.
:type n: int
:param mode: Mode for sampling, defaults to 'random'.
Available modes include: random sampling, 'random';
latin hypercube sampling, 'latin' or 'lh';
chebyshev sampling, 'chebyshev'; grid sampling 'grid'.
:type mode: str, optional.
:param variables: pinn variable to be sampled, defaults to 'all'.
:type variables: str or list[str], optional.
:return: Sample points.
:rtype: list[torch.Tensor]
"""
pairs = [(k, v) for k, v in self.range_.items() if k in variables]
keys, values = map(list, zip(*pairs))
bounds = torch.tensor(values)
result = self._sample_range(n, mode, bounds)
result = result.as_subclass(LabelTensor)
result.labels = keys
for variable in variables:
if variable in self.fixed_.keys():
value = self.fixed_[variable]
pts_variable = torch.tensor([[value]]).repeat(
result.shape[0], 1)
pts_variable = pts_variable.as_subclass(LabelTensor)
pts_variable.labels = [variable]
result = result.append(pts_variable, mode='std')
return result
def _single_points_sample(n, variables):
"""Sample a single point in one dimension.
:param n: Number of points to sample.
:type n: int
:param variables: Variables to sample from.
:type variables: list[str]
:return: Sample points.
:rtype: list[torch.Tensor]
"""
tmp = []
for variable in variables:
if variable in self.fixed_.keys():
value = self.fixed_[variable]
pts_variable = torch.tensor([[value]]).repeat(n, 1)
pts_variable = pts_variable.as_subclass(LabelTensor)
pts_variable.labels = [variable]
tmp.append(pts_variable)
result = tmp[0]
for i in tmp[1:]:
result = result.append(i, mode='std')
return result
if self.fixed_ and (not self.range_):
return _single_points_sample(n, variables)
if variables == 'all':
variables = list(self.range_.keys()) + list(self.fixed_.keys())
if mode in ['grid', 'chebyshev']:
return _1d_sampler(n, mode, variables)
elif mode in ['random', 'lh', 'latin']:
return _Nd_sampler(n, mode, variables)
else:
raise ValueError(f'mode={mode} is not valid.')
def is_inside(self, point, check_border=False):
"""Check if a point is inside the ellipsoid.
: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.
:type check_border: bool
:return: Returning True if the point is inside, False otherwise.
:rtype: bool
"""
is_inside = []
for variable, bound in self.range_.items():
if variable in point.labels:
if bound[0] <= point.extract([variable]) <= bound[1]:
is_inside.append(True)
else:
is_inside.append(False)
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

View File

@@ -0,0 +1,27 @@
"""Module for Location class."""
from .location import Location
from ..label_tensor import LabelTensor
class Difference(Location):
"""
"""
def __init__(self, first, second):
self.first = first
self.second = second
def sample(self, n, mode ='random', variables='all'):
"""
"""
assert mode is '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])
import torch
return LabelTensor(torch.tensor(samples), labels=['x', 'y'])

268
pina/geometry/ellipsoid.py Normal file
View File

@@ -0,0 +1,268 @@
import torch
from .location import Location
from ..label_tensor import LabelTensor
class EllipsoidDomain(Location):
"""PINA implementation of Ellipsoid domain."""
def __init__(self, ellipsoid_dict, sample_surface=False):
"""PINA implementation of Ellipsoid domain.
:param ellipsoid_dict: A dictionary with dict-key a string representing
the input variables for the pinn, and dict-value a list with
the domain extrema.
:type ellipsoid_dict: dict
:param sample_surface: A variable for choosing sample strategies. If
`sample_surface=True` only samples on the ellipsoid surface
frontier are taken. If `sample_surface=False` only samples on
the ellipsoid interior are taken, defaults to False.
:type sample_surface: bool, optional
.. warning::
Sampling for dimensions greater or equal to 10 could result
in a shrinking of the ellipsoid, which degrades the quality
of the samples. For dimensions higher than 10, other algorithms
for sampling should be used, such as: Dezert, Jean, and Christian
Musso. "An efficient method for generating points uniformly
distributed in hyperellipsoids." Proceedings of the Workshop on
Estimation, Tracking and Fusion: A Tribute to Yaakov Bar-Shalom.
Vol. 7. No. 8. 2001.
:Example:
>>> spatial_domain = Ellipsoid({'x':[-1, 1], 'y':[-1,1]})
"""
self.fixed_ = {}
self.range_ = {}
self._centers = None
self._axis = None
if not isinstance(sample_surface, bool):
raise ValueError('sample_surface must be bool type.')
self._sample_surface = sample_surface
for k, v in ellipsoid_dict.items():
if isinstance(v, (int, float)):
self.fixed_[k] = v
elif isinstance(v, (list, tuple)) and len(v) == 2:
self.range_[k] = v
else:
raise TypeError
# perform operation only for not fixed variables (if any)
if self.range_:
# convert dict vals to torch [dim, 2] matrix
list_dict_vals = list(self.range_.values())
tmp = torch.tensor(list_dict_vals, dtype=torch.float)
# get the ellipsoid center
normal_basis = torch.eye(len(list_dict_vals))
centers = torch.diag(normal_basis * tmp.mean(axis=1))
# get the ellipsoid axis
ellipsoid_axis = (tmp - centers.reshape(-1, 1))[:, -1]
# save elipsoid axis and centers as dict
self._centers = dict(zip(self.range_.keys(), centers.tolist()))
self._axis = dict(zip(self.range_.keys(), ellipsoid_axis.tolist()))
@property
def variables(self):
"""Spatial variables.
:return: Spatial variables defined in '__init__()'
:rtype: list[str]
"""
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.
: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.
:type check_border: bool
:return: Returning True if the point is inside, False otherwise.
:rtype: bool
"""
if not isinstance(point, LabelTensor):
raise ValueError('point expected to be LabelTensor.')
# get axis ellipse
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()))
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_sq.labels = point.labels
# calculate ellispoid equation
eqn = torch.sum(point_sq.extract(ax_sq.labels) / ax_sq) - 1.
if check_border:
return bool(eqn <= 0)
return bool(eqn < 0)
def _sample_range(self, n, mode, variables):
"""Rescale the samples to the correct bounds.
:param n: Number of points to sample in the ellipsoid.
:type n: int
:param mode: Mode for sampling, defaults to 'random'.
Available modes include: random sampling, 'random'.
:type mode: str, optional
:param variables: Variables to be rescaled in the samples.
:type variables: torch.Tensor
:return: Rescaled sample points.
:rtype: torch.Tensor
"""
# =============== For Developers ================ #
#
# The sampling startegy used is fairly simple.
# For all `mode`s first we sample from the unit
# sphere and then we scale and shift according
# to self._axis.values() and self._centers.values().
#
# =============================================== #
# get dimension
dim = len(variables)
# get values center
pairs_center = [(k, v) for k, v in self._centers.items()
if k in variables]
_, values_center = map(list, zip(*pairs_center))
values_center = torch.tensor(values_center)
# get values axis
pairs_axis = [(k, v) for k, v in self._axis.items() if k in variables]
_, values_axis = map(list, zip(*pairs_axis))
values_axis = torch.tensor(values_axis)
# Sample in the unit sphere
if mode == 'random':
# 1. Sample n points from the surface of a unit sphere
# 2. Scale each dimension using torch.rand()
# (a random number between 0-1) so that it lies within
# the sphere, only if self._sample_surface=False
# 3. Multiply with self._axis.values() to make it ellipsoid
# 4. Shift the mean of the ellipse by adding self._centers.values()
# step 1.
pts = torch.randn(size=(n, dim))
pts = pts / torch.linalg.norm(pts, axis=-1).view((n, 1))
if not self._sample_surface: # step 2.
scale = torch.rand((n, 1))
pts = pts * scale
# step 3. and 4.
pts *= values_axis
pts += values_center
return pts
def sample(self, n, mode='random', variables='all'):
"""Sample routine.
:param n: Number of points to sample in the ellipsoid.
:type n: int
:param mode: Mode for sampling, defaults to 'random'.
Available modes include: random sampling, 'random'.
:type mode: str, optional
:param variables: pinn variable to be sampled, defaults to 'all'.
:type variables: str or list[str], optional
:Example:
>>> elips = Ellipsoid({'x':[1, 0], 'y':1})
>>> elips.sample(n=6)
tensor([[0.4872, 1.0000],
[0.2977, 1.0000],
[0.0422, 1.0000],
[0.6431, 1.0000],
[0.7272, 1.0000],
[0.8326, 1.0000]])
"""
def _Nd_sampler(n, mode, variables):
"""Sample all the variables together
:param n: Number of points to sample.
:type n: int
:param mode: Mode for sampling, defaults to 'random'.
Available modes include: random sampling, 'random';
latin hypercube sampling, 'latin' or 'lh';
chebyshev sampling, 'chebyshev'; grid sampling 'grid'.
:type mode: str, optional.
:param variables: pinn variable to be sampled, defaults to 'all'.
:type variables: str or list[str], optional.
:return: Sample points.
:rtype: list[torch.Tensor]
"""
pairs = [(k, v) for k, v in self.range_.items() if k in variables]
keys, _ = map(list, zip(*pairs))
result = self._sample_range(n, mode, keys)
result = result.as_subclass(LabelTensor)
result.labels = keys
for variable in variables:
if variable in self.fixed_.keys():
value = self.fixed_[variable]
pts_variable = torch.tensor([[value]
]).repeat(result.shape[0], 1)
pts_variable = pts_variable.as_subclass(LabelTensor)
pts_variable.labels = [variable]
result = result.append(pts_variable, mode='std')
return result
def _single_points_sample(n, variables):
"""Sample a single point in one dimension.
:param n: Number of points to sample.
:type n: int
:param variables: Variables to sample from.
:type variables: list[str]
:return: Sample points.
:rtype: list[torch.Tensor]
"""
tmp = []
for variable in variables:
if variable in self.fixed_.keys():
value = self.fixed_[variable]
pts_variable = torch.tensor([[value]]).repeat(n, 1)
pts_variable = pts_variable.as_subclass(LabelTensor)
pts_variable.labels = [variable]
tmp.append(pts_variable)
result = tmp[0]
for i in tmp[1:]:
result = result.append(i, mode='std')
return result
if self.fixed_ and (not self.range_):
return _single_points_sample(n, variables)
if variables == 'all':
variables = list(self.range_.keys()) + list(self.fixed_.keys())
if mode in ['random']:
return _Nd_sampler(n, mode, variables)
else:
raise ValueError(f'mode={mode} is not valid.')

17
pina/geometry/location.py Normal file
View File

@@ -0,0 +1,17 @@
"""Module for Location class."""
from abc import ABCMeta, abstractmethod
class Location(metaclass=ABCMeta):
"""
Abstract Location class.
Any geometry entity should inherit from this class.
"""
@abstractmethod
def sample(self):
"""
Abstract method for sampling a point from the location. To be
implemented in the child class.
"""
pass