Commit a11decc5 authored by dboe's avatar dboe
Browse files

Merge branch 'dev'

parents 073c0b7e 269b6256
Pipeline #108029 passed with stages
in 53 seconds
import numpy as np
from tempfile import NamedTemporaryFile
# pylint:disable=missing-function-docstring,missing-module-docstring,missing-class-docstring,invalid-name,protected-access
import pickle
import pathlib
import unittest
import uuid
from tempfile import NamedTemporaryFile
from copy import deepcopy
import numpy as np
import tfields
ATOL = 1e-8
class Base_Check(object):
# pylint:disable=no-member
class Base_Check:
def demand_equal(self, other):
self.assertIsInstance(other, type(self._inst))
......@@ -26,8 +31,6 @@ class Base_Check(object):
self.demand_equal(reloaded)
def test_deep_copy(self):
from copy import deepcopy
other = deepcopy(self._inst)
self.demand_deep_copy(other)
......@@ -61,6 +64,7 @@ class AbstractNdarray_Check(Base_Check):
pass
# pylint:disable=no-member
class Tensors_Check(AbstractNdarray_Check):
"""
Testing derivatives of Tensors
......@@ -100,6 +104,7 @@ class Tensors_Check(AbstractNdarray_Check):
tensors = np.array(self._inst)
if len(self._inst) > 0:
# pylint:disable=unsubscriptable-object
item = self._inst[index]
self.assertTrue(np.array_equal(item, tensors[index]))
self.assertIsInstance(item, check_type)
......@@ -165,6 +170,7 @@ class Tensors_Check(AbstractNdarray_Check):
self.assertTrue(value)
# pylint:disable=no-member
class TensorFields_Check(Tensors_Check):
def test_fields(self):
self.assertIsNotNone(self._inst.fields)
......@@ -173,15 +179,21 @@ class TensorFields_Check(Tensors_Check):
self.assertTrue(isinstance(self._inst.fields, list))
self.assertTrue(len(self._inst.fields) == len(self._fields))
for field, target_field in zip(self._inst.fields, self._fields):
self.assertTrue(np.array_equal(field, target_field))
# fields are copied not reffered by a pointer
self.assertFalse(field is target_field)
self.check_fields_equal(self._inst.fields, self._fields)
def check_fields_equal(self, fields_a, fields_b):
for field, target_field in zip(fields_a, fields_b):
self.assertTrue(np.array_equal(field, target_field))
# The below test is purposely disabled. It is better as a philosphy to not leave
# setitem values mutable.
# # fields are copied not reffered by a pointer
# self.assertFalse(field is target_field)
def demand_index_equal(self, index, check_type):
super().demand_index_equal(index, check_type)
if len(self._inst) > 0:
# pylint:disable=unsubscriptable-object
item = self._inst[index]
for i, field in enumerate(self._inst.fields):
if check_type == "type":
......@@ -199,6 +211,44 @@ class TensorFields_Check(Tensors_Check):
for i in range(len(self._inst.fields)):
self.assertIsNot(self._inst.fields[i], other.fields[i])
def test_list_like_field(self):
if self._inst.fields:
fields = self._inst.fields
self._inst.fields = []
for i, field in enumerate(fields):
self._inst.fields.append(field)
# indexing
self.assertTrue(self._inst.fields[i].equal(field))
self.check_fields_equal(fields, self._inst.fields)
def test_field_name_getitem(self):
if self._inst.fields:
for field in self._inst.fields:
if field.name is not None:
self.assertTrue(self._inst.fields[field.name].equal(field))
def test_field_name_setitem(self):
if self._inst.fields:
fields = self._inst.fields
self._inst.fields = []
for i, field in enumerate(fields):
if field.name is not None:
name = field.name
else:
name = str(uuid.uuid4())
field.name = name
# setitem via fields
self._inst.fields[name] = field
field_item = self._inst.fields[name]
self.assertTrue(field_item.equal(field))
self.assertTrue(self._inst.fields[i].equal(field_item))
self.assertEqual(field_item.name, name)
self.check_fields_equal(fields, self._inst.fields)
class TensorMaps_Check(TensorFields_Check):
def test_maps(self):
......@@ -221,9 +271,9 @@ class TensorMaps_Check(TensorFields_Check):
self.assertIsNot(self._inst.maps[i], other.maps[i])
"""
EMPTY TESTS
"""
#############
# EMPTY TESTS
#############
class Tensors_Empty_Test(Tensors_Check, unittest.TestCase):
......@@ -245,7 +295,7 @@ class TensorMaps_Empty_Test(TensorMaps_Check, unittest.TestCase):
self._maps_fields = []
class TensorFields_Test(Tensors_Check, unittest.TestCase):
class TensorFields_Test(TensorFields_Check, unittest.TestCase):
def setUp(self):
base = [(-5, 5, 11)] * 3
self._fields = [
......@@ -259,7 +309,7 @@ class TensorFields_Test(Tensors_Check, unittest.TestCase):
self.assertTrue(self._fields[1].coord_sys, "cartesian")
class TensorMaps_Test(Tensors_Check, unittest.TestCase):
class TensorMaps_Test(TensorFields_Check, unittest.TestCase):
def setUp(self):
base = [(-1, 1, 3)] * 3
tensors = tfields.Tensors.grid(*base)
......
# pylint:disable=missing-function-docstring,missing-module-docstring,missing-class-docstring,invalid-name
import unittest
import itertools
import numpy as np
......@@ -43,7 +44,7 @@ class TensorGrid_Test(unittest.TestCase, TensorFields_Check):
def check_filled(self, tg):
self.assertTrue(tg.equal(self.res))
self.assertListEqual(tg.base_vectors, self.bv)
self.assertListEqual(list(tg.base_num_tuples()), self.bv)
def check_empty(self, tg):
self.assertTrue(tg.is_empty())
......@@ -62,9 +63,17 @@ class TensorGrid_Test(unittest.TestCase, TensorFields_Check):
*self.bv, iter_order=self.iter_order, coord_sys=tfields.bases.CYLINDER
)
self.assertEqual(tg.coord_sys, tfields.bases.CYLINDER)
self.assertEqual(tg.base_vectors.coord_sys, tfields.bases.CYLINDER)
tge = tg.explicit()
self.assertEqual(tge.coord_sys, tfields.bases.CYLINDER)
# transformation does not change the base_vector coord_sys
tg.transform(tfields.bases.CARTESIAN)
self.assertEqual(tg.coord_sys, tfields.bases.CARTESIAN)
self.assertEqual(tg.base_vectors.coord_sys, tfields.bases.CYLINDER)
tge = tg.explicit()
self.assertEqual(tge.coord_sys, tfields.bases.CARTESIAN)
def test_getitem(self):
tg = TensorGrid.empty(*self.bv, iter_order=self.iter_order)
self.check_filled(tg[:])
......@@ -111,15 +120,29 @@ 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()
t2 = TensorGrid_Test_Permutation1()
t2.setUp()
# pylint:disable=protected-access
tg1 = t1._inst
tg1.fields.append(tfields.Tensors(tg1))
tg1_original = tg1.copy()
# pylint:disable=protected-access
tg2 = t2._inst
tg2.fields.append(tfields.Tensors(tg2))
......@@ -137,15 +160,12 @@ class TensorGrid_Test_IO_Change(unittest.TestCase):
def test_lib_multi(self):
bv_lengths = [2, 3, 4]
# pylint:disable=unused-variable
for iter_order in itertools.permutations([0, 1, 2]):
# pylint:disable=unused-variable
for new_iter_order in itertools.permutations([0, 1, 2]):
forw = grid.change_iter_order(bv_lengths, [0, 1, 2], [2, 0, 1])
back = grid.change_iter_order(bv_lengths, [2, 0, 1], [0, 1, 2])
self.assertTrue(
(np.array(forw)[back] == np.arange(np.prod(bv_lengths))).all()
)
if __name__ == "__main__":
TensorGrid_Test().test_tensor_grid()
# unittest.main()
......@@ -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"