Commit 269b6256 authored by dboe's avatar dboe
Browse files

added PHYSICAL_CYLINDER, transform_field and lib.io

parent 6b158fd3
Pipeline #108026 failed with stage
in 25 seconds
......@@ -120,6 +120,18 @@ class TensorGrid_Test_Permutation1(TensorGrid_Test):
class TensorGrid_Test_IO_Change(unittest.TestCase):
def test_copy_constructor(self):
tgt = TensorGrid_Test()
tgt.setUp()
tg_orig = tgt._inst
tg_copy = tfields.TensorGrid(tg_orig)
tgt.check_filled(tg_copy)
tg_copy_new_num = tfields.TensorGrid(tg_orig, num=[1, 1, 1])
self.assertListEqual(list(tg_copy_new_num.num), [1, 1, 1])
self.assertEqual(len(tg_copy_new_num), 1)
def test_change_iter_order(self):
t1 = TensorGrid_Test()
t1.setUp()
......
......@@ -7,7 +7,24 @@ Mail: daniel.boeckenhoff@ipp.mpg.de
part of tfields library
Tools for sympy coordinate transformation
"""
from tfields.bases.bases import get_coord_system, get_coord_system_name, lambdified_trafo, transform # NOQA
from tfields.bases.bases import (
get_coord_system,
get_coord_system_name,
lambdified_trafo,
transform,
) # NOQA
from tfields.bases import manifold_3 # NOQA
from tfields.bases.manifold_3 import CARTESIAN, CYLINDER, SPHERICAL # NOQA
from tfields.bases.manifold_3 import cartesian, cylinder, spherical # NOQA
from tfields.bases.manifold_3 import (
CARTESIAN,
CYLINDER,
SPHERICAL,
NATURAL_CARTESIAN,
PHYSICAL_CYLINDER,
) # NOQA
from tfields.bases.manifold_3 import (
cartesian,
cylinder,
spherical,
natural_cartesian,
physical_cylinder,
) # NOQA
......@@ -22,15 +22,17 @@ def get_coord_system(base):
Return:
sympy.diffgeom.get_coord_system
"""
if isinstance(base, string_types) or (isinstance(base, np.ndarray)
and base.dtype.kind in {'U', 'S'}):
if isinstance(base, string_types) or (
isinstance(base, np.ndarray) and base.dtype.kind in {"U", "S"}
):
base = getattr(tfields.bases, str(base))
if not isinstance(base, sympy.diffgeom.CoordSystem):
bse_tpe = type(base)
expctd_tpe = type(sympy.diffgeom.CoordSystem)
raise TypeError("Wrong type of coord_system base {bse_tpe}. "
"Expected {expctd_tpe}"
.format(**locals()))
raise TypeError(
"Wrong type of coord_system base {bse_tpe}. "
"Expected {expctd_tpe}".format(**locals())
)
return base
......@@ -42,7 +44,7 @@ def get_coord_system_name(base):
str: name of base
"""
if isinstance(base, sympy.diffgeom.CoordSystem):
base = getattr(base, 'name')
base = getattr(base, "name")
# if not (isinstance(base, string_types) or base is None):
# baseType = type(base)
# raise ValueError("Coordinate system must be string_type."
......@@ -76,14 +78,15 @@ def lambdified_trafo(base_old, base_new):
"""
coords = tuple(base_old.coord_function(i) for i in range(base_old.dim))
f = sympy.lambdify(coords,
base_old.coord_tuple_transform_to(base_new,
list(coords)),
modules='numpy')
f = sympy.lambdify(
coords,
base_old.coord_tuple_transform_to(base_new, list(coords)),
modules="numpy",
)
return f
def transform(array, base_old, base_new):
def transform(array, base_old, base_new, **kwargs):
"""
Transform the input array in place
Args:
......@@ -119,31 +122,35 @@ def transform(array, base_old, base_new):
"""
base_old = get_coord_system(base_old)
base_new = get_coord_system(base_new)
if base_new not in base_old.transforms:
for baseTmp in base_new.transforms:
if baseTmp in base_old.transforms:
transform(array, base_old, baseTmp)
transform(array, baseTmp, base_new)
return
raise ValueError("Trafo not found.")
# very fast trafos in numpy only
short_trafo = None
try:
short_trafo = getattr(base_old, 'to_{base_new.name}'.format(**locals()))
short_trafo = getattr(base_old, f"to_{base_new.name}".format(**locals()))
except AttributeError:
pass
if short_trafo:
short_trafo(array)
short_trafo(array, **kwargs)
return
if base_new not in base_old.transforms:
for baseTmp in base_new.transforms:
if baseTmp in base_old.transforms:
transform(array, base_old, baseTmp, **kwargs)
transform(array, baseTmp, base_new, **kwargs)
return
raise ValueError(f"Transformation {base_old} -> {base_new} not found.")
# trafo via lambdified sympy expressions
trafo = tfields.bases.lambdified_trafo(base_old, base_new)
with warnings.catch_warnings():
warnings.filterwarnings('ignore', message="invalid value encountered in double_scalars")
warnings.filterwarnings(
"ignore", message="invalid value encountered in double_scalars"
)
array[:] = np.concatenate([trafo(*coords).T for coords in array])
if __name__ == '__main__': # pragma: no cover
if __name__ == "__main__": # pragma: no cover
import doctest
doctest.testmod()
......@@ -4,80 +4,96 @@ import sympy.diffgeom
import numpy as np
import warnings
CARTESIAN = 'cartesian'
CYLINDER = 'cylinder'
SPHERICAL = 'spherical'
m = sympy.diffgeom.Manifold('M', 3)
patch = sympy.diffgeom.Patch('P', m)
# cartesian
x, y, z = sympy.symbols('x, y, z')
cartesian = sympy.diffgeom.CoordSystem(CARTESIAN, patch, ['x', 'y', 'z'])
# cylinder
R, Phi, Z = sympy.symbols('R, Phi, Z')
cylinder = sympy.diffgeom.CoordSystem(CYLINDER, patch, ['R', 'Phi', 'Z'])
cylinder.connect_to(cartesian,
[R, Phi, Z],
[R * sympy.cos(Phi),
R * sympy.sin(Phi),
Z],
inverse=False,
fill_in_gaps=False)
CARTESIAN = "cartesian"
CYLINDER = "cylinder"
SPHERICAL = "spherical"
m = sympy.diffgeom.Manifold("M", 3)
patch = sympy.diffgeom.Patch("P", m)
#############
# cartesian #
#############
x, y, z = sympy.symbols("x, y, z")
cartesian = sympy.diffgeom.CoordSystem(CARTESIAN, patch, ["x", "y", "z"])
############
# cylinder #
############
R, Phi, Z = sympy.symbols("R, Phi, Z")
cylinder = sympy.diffgeom.CoordSystem(CYLINDER, patch, ["R", "Phi", "Z"])
cylinder.connect_to(
cartesian,
[R, Phi, Z],
[R * sympy.cos(Phi), R * sympy.sin(Phi), Z],
inverse=False,
fill_in_gaps=False,
)
def cylinder_to_cartesian(array):
rPoints = np.copy(array[:, 0])
phiPoints = np.copy(array[:, 1])
array[:, 0] = rPoints * np.cos(phiPoints)
array[:, 1] = rPoints * np.sin(phiPoints)
r_points = np.copy(array[:, 0])
phi_points = np.copy(array[:, 1])
array[:, 0] = r_points * np.cos(phi_points)
array[:, 1] = r_points * np.sin(phi_points)
cylinder.to_cartesian = cylinder_to_cartesian
cartesian.connect_to(cylinder,
[x, y, z],
[sympy.sqrt(x**2 + y**2),
sympy.Piecewise(
(sympy.pi, ((x < 0) & sympy.Eq(y, 0))),
(0., sympy.Eq(sympy.sqrt(x**2 + y**2), 0)),
(sympy.sign(y) * sympy.acos(x / sympy.sqrt(x**2 + y**2)), True)),
z],
inverse=False,
fill_in_gaps=False)
cartesian.connect_to(
cylinder,
[x, y, z],
[
sympy.sqrt(x ** 2 + y ** 2),
sympy.Piecewise(
(sympy.pi, ((x < 0) & sympy.Eq(y, 0))),
(0.0, sympy.Eq(sympy.sqrt(x ** 2 + y ** 2), 0)),
(sympy.sign(y) * sympy.acos(x / sympy.sqrt(x ** 2 + y ** 2)), True),
),
z,
],
inverse=False,
fill_in_gaps=False,
)
def cartesian_to_cylinder(array):
x_vals = np.copy(array[:, 0])
y_vals = np.copy(array[:, 1])
problem_phi_indices = np.where((x_vals < 0) & (y_vals == 0))
array[:, 0] = np.sqrt(x_vals**2 + y_vals**2) # r
np.seterr(divide='ignore', invalid='ignore')
array[:, 0] = np.sqrt(x_vals ** 2 + y_vals ** 2) # r
np.seterr(divide="ignore", invalid="ignore")
# phi
array[:, 1] = np.sign(y_vals) * np.arccos(x_vals / array[:, 0])
np.seterr(divide='warn', invalid='warn')
np.seterr(divide="warn", invalid="warn")
array[:, 1][problem_phi_indices] = np.pi
tfields.lib.util.convert_nan(array, 0.) # for cases like cartesian 0, 0, 1
tfields.lib.util.convert_nan(array, 0.0) # for cases like cartesian 0, 0, 1
cartesian.to_cylinder = cartesian_to_cylinder
# spherical
r, phi, theta = sympy.symbols('r, phi, theta')
spherical = sympy.diffgeom.CoordSystem(SPHERICAL, patch, ['r', 'phi', 'theta'])
spherical.connect_to(cartesian,
[r, phi, theta],
[r * sympy.cos(phi - np.pi) * sympy.sin(theta - np.pi / 2),
r * sympy.sin(phi - np.pi) * sympy.sin(theta - np.pi / 2),
r * sympy.cos(theta - np.pi / 2)],
# [r * sympy.cos(phi) * sympy.sin(theta),
# r * sympy.sin(phi) * sympy.sin(theta),
# r * sympy.cos(theta)],
inverse=False,
fill_in_gaps=False)
#############
# spherical #
#############
r, phi, theta = sympy.symbols("r, phi, theta")
spherical = sympy.diffgeom.CoordSystem(SPHERICAL, patch, ["r", "phi", "theta"])
spherical.connect_to(
cartesian,
[r, phi, theta],
[
r * sympy.cos(phi - np.pi) * sympy.sin(theta - np.pi / 2),
r * sympy.sin(phi - np.pi) * sympy.sin(theta - np.pi / 2),
r * sympy.cos(theta - np.pi / 2),
],
# [r * sympy.cos(phi) * sympy.sin(theta),
# r * sympy.sin(phi) * sympy.sin(theta),
# r * sympy.cos(theta)],
inverse=False,
fill_in_gaps=False,
)
def spherical_to_cartesian(array):
......@@ -97,26 +113,31 @@ def spherical_to_cartesian(array):
spherical.to_cartesian = spherical_to_cartesian
cartesian.connect_to(spherical,
[x, y, z],
[sympy.sqrt(x**2 + y**2 + z**2),
sympy.Piecewise(
(sympy.pi, ((x < 0) & sympy.Eq(y, 0))),
(0., sympy.Eq(sympy.sqrt(x**2 + y**2), 0)),
(sympy.sign(y) * sympy.acos(x / sympy.sqrt(x**2 + y**2)), True)),
sympy.Piecewise(
(0., sympy.Eq(x**2 + y**2 + z**2, 0)),
(sympy.atan2(z, sympy.sqrt(x**2 + y**2)), True)),
# sympy.Piecewise(
# (0., (sympy.Eq(x, 0) & sympy.Eq(y, 0))),
# (sympy.atan(y / x), True)),
# sympy.Piecewise(
# (0., sympy.Eq(x**2 + y**2 + z**2, 0)),
# # (0., (sympy.Eq(x, 0) & sympy.Eq(y, 0) & sympy.Eq(z, 0))),
# (sympy.acos(z / sympy.sqrt(x**2 + y**2 + z**2)), True))
],
inverse=False,
fill_in_gaps=False)
cartesian.connect_to(
spherical,
[x, y, z],
[
sympy.sqrt(x ** 2 + y ** 2 + z ** 2),
sympy.Piecewise(
(sympy.pi, ((x < 0) & sympy.Eq(y, 0))),
(0.0, sympy.Eq(sympy.sqrt(x ** 2 + y ** 2), 0)),
(sympy.sign(y) * sympy.acos(x / sympy.sqrt(x ** 2 + y ** 2)), True),
),
sympy.Piecewise(
(0.0, sympy.Eq(x ** 2 + y ** 2 + z ** 2, 0)),
(sympy.atan2(z, sympy.sqrt(x ** 2 + y ** 2)), True),
),
# sympy.Piecewise(
# (0., (sympy.Eq(x, 0) & sympy.Eq(y, 0))),
# (sympy.atan(y / x), True)),
# sympy.Piecewise(
# (0., sympy.Eq(x**2 + y**2 + z**2, 0)),
# # (0., (sympy.Eq(x, 0) & sympy.Eq(y, 0) & sympy.Eq(z, 0))),
# (sympy.acos(z / sympy.sqrt(x**2 + y**2 + z**2)), True))
],
inverse=False,
fill_in_gaps=False,
)
def cartesian_to_spherical(array):
......@@ -126,25 +147,25 @@ def cartesian_to_spherical(array):
phi in [-pi, pi]
theta element [0, pi]: elevation angle defined from - Z-axis up
"""
xy = array[:, 0]**2 + array[:, 1]**2
xy = array[:, 0] ** 2 + array[:, 1] ** 2
# phi for phi between -pi, pi
problemPhiIndices = np.where((array[:, 0] < 0) & (array[:, 1] == 0))
with warnings.catch_warnings():
# python2.7
warnings.filterwarnings('ignore',
message="invalid value encountered in divide")
warnings.filterwarnings("ignore", message="invalid value encountered in divide")
# python3.x
warnings.filterwarnings('ignore',
message="invalid value encountered in true_divide")
warnings.filterwarnings(
"ignore", message="invalid value encountered in true_divide"
)
array[:, 1] = np.sign(array[:, 1]) * np.arccos(array[:, 0] / np.sqrt(xy))
array[:, 1][problemPhiIndices] = np.pi
tfields.lib.util.convert_nan(array, 0.) # for cases like cartesian 0, 0, 1
tfields.lib.util.convert_nan(array, 0.0) # for cases like cartesian 0, 0, 1
# array[:,1] = np.arctan2(array[:,1], array[:,0]) # phi for phi between 0, 2pi
# r
array[:, 0] = np.sqrt(xy + array[:, 2]**2)
array[:, 0] = np.sqrt(xy + array[:, 2] ** 2)
# theta
# for elevation angle defined from - Z-axis up:
......@@ -154,3 +175,36 @@ def cartesian_to_spherical(array):
cartesian.to_spherical = cartesian_to_spherical
#################################################################################
# TODO: the part below is a temporary hack. I do not know yet the maths properly.
# will have to learn that and generalize.
#################################################################################
PHYSICAL_CYLINDER = "physical_cylinder"
NATURAL_CARTESIAN = "natural_cartesian"
natural_cartesian = sympy.diffgeom.CoordSystem(
NATURAL_CARTESIAN, patch, ["x", "y", "z"]
)
physical_cylinder = sympy.diffgeom.CoordSystem(
PHYSICAL_CYLINDER, patch, ["R", "Phi", "Z"]
)
def physical_cylinder_to_natural_cartesian(vector, position=None):
"""
We need to project from PHYSICAL cylinder coordinates to NATURAL cartesian coordinates
"""
with position.tmp_transform(CYLINDER):
sin_phi = np.copy(np.sin(position[:, 1]))
cos_phi = np.copy(np.cos(position[:, 1]))
vec_r = np.copy(vector[:, 0])
vec_phi = np.copy(vector[:, 1])
vector[:, 0] = vec_r * cos_phi - vec_phi * sin_phi
vector[:, 1] = vec_phi * cos_phi - vec_r * sin_phi
physical_cylinder.to_natural_cartesian = physical_cylinder_to_natural_cartesian
......@@ -17,6 +17,7 @@ TODO:
"""
# builtin
import warnings
import inspect
from contextlib import contextmanager
from collections import Counter
from copy import deepcopy
......@@ -287,7 +288,7 @@ class AbstractNdarray(np.ndarray, AbstractObject):
defaults for the attributes in __slots__ will be None
other values will be treaded as defaults to the corresponding
arg at the same position in the __slots__ list.
__slot_dtype__ (List(dtypes)): for the conversion of the
__slot_dtypes__ (List(dtypes)): for the conversion of the
args in __slots__ to numpy arrays. None values mean no
conversion.
__slot_setters__ (List(callable)): Because __slots__ and properties are
......@@ -704,6 +705,89 @@ class Tensors(AbstractNdarray): # pylint: disable=too-many-public-methods
for index in range(len(self)):
yield super(Tensors, self).__getitem__(index).view(Tensors)
@classmethod
def _load_txt(cls, path, set_header: bool = False, **load_kwargs):
"""
Factory method
Given a path to a txt file, construct the object
Args:
path: see :meth:`Tensors.load`
set_header: if Treu, set the header list to the `name` attribute
**loadkwargs: see :meth:`np.loadtxt`
"""
load_txt_keys = [
k
for k, v in inspect.signature(np.loadtxt).parameters.items()
if v.default is not inspect._empty # pylint:disable=protected-access
]
load_txt_kwargs = {}
for key in load_txt_keys:
if key in load_kwargs:
load_txt_kwargs[key] = load_kwargs.pop(key)
# use the same defaults as np.savetxt
load_txt_kwargs.setdefault("comments", "# ") # default is "#"
comments = load_txt_kwargs.get("comments")
load_txt_kwargs.setdefault("delimiter", " ")
skiprows = 0
header = []
with open(rna.path.resolve(path)) as file_:
while True:
line = file_.readline()
if line.startswith(comments):
header.append(line.lstrip(comments).rstrip("\n"))
else:
file_.seek(0)
break
skiprows += 1
load_txt_kwargs["skiprows"] = max(skiprows, load_txt_kwargs.pop("skiprows", 0))
arr = np.loadtxt(path, **load_txt_kwargs)
obj = cls(arr, **load_kwargs)
i = 0
for key, line in zip(cls._iter_slots(), header):
slot, rest = line.split(": ")
tpe, val_str = rest.split(" = ")
if tpe == "NoneType":
val = None
else:
tpe = __builtins__[tpe]
val = tpe(val_str)
setattr(obj, slot, val)
i += 1
header = header[i:]
if set_header:
obj.name = (
obj.name,
header,
) # pylint:disable=attribute-defined-outside-init
return obj
def _save_txt(self, path, **kwargs):
"""
Save as text file.
Args:
**kwargs passed to np.savetxt.
"""
header = kwargs.get("header", [])
if isinstance(header, dict):
header = [
f"{key}: {type(value).__name__} = {value}"
for key, value in header.items()
]
if isinstance(header, list):
# statictyping like attribute saving
header = [
f"{key}: {type(getattr(self, key)).__name__} = {getattr(self, key)}"
for key in self._iter_slots()
] + header
kwargs["header"] = "\n".join(header)
np.savetxt(path, self, **kwargs)
@classmethod
def merged(cls, *objects, **kwargs):
"""
......@@ -925,7 +1009,7 @@ class Tensors(AbstractNdarray): # pylint: disable=too-many-public-methods
"""
return dim(self)
def transform(self, coord_sys):
def transform(self, coord_sys, **kwargs):
"""
Args:
coord_sys (str)
......@@ -1000,7 +1084,7 @@ class Tensors(AbstractNdarray): # pylint: disable=too-many-public-methods
# already correct
return
tfields.bases.transform(self, self.coord_sys, coord_sys)
tfields.bases.transform(self, self.coord_sys, coord_sys, **kwargs)
# self[:] = tfields.bases.transform(self, self.coord_sys, coord_sys)
self.coord_sys = coord_sys # pylint: disable=attribute-defined-outside-init
......@@ -1695,6 +1779,8 @@ class Tensors(AbstractNdarray): # pylint: disable=too-many-public-methods
def plot(self, *args, **kwargs):
"""
Generic plotting method of Tensors.
Forwarding to rna.plotting.plot_tensor
"""
artist = rna.plotting.plot_tensor(
......@@ -1967,6 +2053,15 @@ class TensorFields(Tensors):
else:
return inst
def transform_field(self, coord_sys, field_index=0):
"""
Transform the field to the coordinate system of choice.
NOTE: This is not yet any way generic!!! Have a look at Einsteinpy and actual status of
sympy for further implementation
"""
self.fields[field_index].transform(coord_sys, position=self)
@property
def names(self):
"""
......@@ -2002,8 +2097,7 @@ class TensorFields(Tensors):
Args:
other (iterable)
optional:
see Tensors.equal
**kwargs: see Tensors.equal
"""
if not issubclass(type(other), Tensors):
return super(TensorFields, self).equal(other, **kwargs)
......@@ -2032,11 +2126,11 @@ class TensorFields(Tensors):
def plot(self, *args, **kwargs):
"""
Plotting the tensor field.
Generic plotting method of TensorFields.
Args:
field_index: index of the field to plot (as quiver by default)
normalize: If true, normalize the field vectors to show only the direction
normalize: If True, normalize the field vectors to show only the direction
color: additional str argument 'norm' added. If color="norm", color with the norm.
"""
field_index = kwargs.pop("field_index", None)
......@@ -2053,13 +2147,14 @@ class TensorFields(Tensors):
color = kwargs.get("color", None)
field = self.fields[field_index].copy()
if self.dim == field.dim:
field.transform(self.coord_sys)
else:
logging.debug(
"Careful: Plotting tensors with field of"
"different dimension. No coord_sys check performed."
)
# this tranfomration is compmlicated. Transforming the field is not yet understood
# if self.dim == field.dim:
# field.transform(self.coord_sys)
# else:
# logging.debug(
# "Careful: Plotting tensors with field of"
# "different dimension. No coord_sys check performed."
# )
if color == "norm":
norm = field.norm()
...