Commit cd583406 authored by Daniel Boeckenhoff's avatar Daniel Boeckenhoff
Browse files

Merge branch 'master' of gitlab.mpcdf.mpg.de:dboe/tfields

parents b460bc66 1eff624f
# Byte-compiled / optimized / DLL files
__pycache__/
*.py[cod]
*$py.class
# C extensions
*.so
# Distribution / packaging
.Python
env/
build/
develop-eggs/
dist/
downloads/
eggs/
.eggs/
lib/
lib64/
parts/
sdist/
var/
*.egg-info/
.installed.cfg
*.egg
# PyInstaller
# Usually these files are written by a python script from a template
# before PyInstaller builds the exe, so as to inject date/other infos into it.
*.manifest
*.spec
# Installer logs
pip-log.txt
pip-delete-this-directory.txt
# Unit test / coverage reports
htmlcov/
.tox/
.coverage
.coverage.*
.cache
nosetests.xml
coverage.xml
*,cover
.hypothesis/
# pyc files:
*.pyc
# Translations
*.mo
*.pot
# Django stuff:
*.log
local_settings.py
# Flask instance folder
instance/
# Scrapy stuff:
.scrapy
# Sphinx documentation
docs/_build/
# PyBuilder
target/
# IPython Notebook
.ipynb_checkpoints
# pyenv
.python-version
# celery beat schedule file
celerybeat-schedule
# dotenv
.env
# virtualenv
venv/
ENV/
# Spyder project settings
.spyderproject
# vim files
*.sw*
# folder
tmp/
.idea
from . import core
from . import bases
from . import lib
from .lib import *
# __all__ = ['core', 'points3D']
from .core import Tensors, TensorFields, TensorMaps
# from .points3D import Points3D
from .points3D import Points3D
from .mask import getMask
# methods:
from .mask import getMask # NOQA
from .lib import * # NOQA
# classes:
from .points3D import Points3D # NOQA
from .mesh3D import Mesh3D # NOQA
from .triangles3D import Triangles3D # NOQA
from .scalarField3D import ScalarField3D # NOQA
from .vectorField3D import VectorField3D # NOQA
from .planes3D import Planes3D # NOQA
......@@ -4,8 +4,9 @@ Run doctests on module call
def runDoctests():
import doctest
import pathlib
for f in list(pathlib.Path(__file__).parent.glob('**/*.py')):
doctest.testfile(f.name) # , verbose=True, optionflags=doctest.ELLIPSIS)
parent = pathlib.Path(__file__).parent
for f in list(parent.glob('**/*.py')):
doctest.testfile(str(f.relative_to(parent))) # , verbose=True, optionflags=doctest.ELLIPSIS)
if __name__ == '__main__':
......
......@@ -14,68 +14,17 @@ import sympy.diffgeom
from six import string_types
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.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.acos(x / sympy.sqrt(x**2 + y**2)), True)),
z],
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) * sympy.sin(theta),
r * sympy.sin(phi) * sympy.sin(theta),
r * sympy.cos(theta)],
inverse=False,
fill_in_gaps=False)
cartesian.connect_to(spherical,
[x, y, z],
[sympy.sqrt(x**2 + y**2 + z**2),
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 getCoordSystem(base):
from . import manifold_3
from .manifold_3 import CARTESIAN, CYLINDER, SPHERICAL
from .manifold_3 import cartesian, cylinder, spherical
def get_coord_system(base):
"""
Args:
base (str or sympy.diffgeom.getCoordSystem)
base (str or sympy.diffgeom.get_coord_system)
Return:
sympy.diffgeom.getCoordSystem
sympy.diffgeom.get_coord_system
"""
if isinstance(base, string_types):
base = getattr(tfields.bases, base)
......@@ -84,10 +33,10 @@ def getCoordSystem(base):
return base
def getCoordSystemName(base):
def get_coord_system_name(base):
"""
Args:
base (str or sympy.diffgeom.getCoordSystem)
base (str or sympy.diffgeom.get_coord_system)
Returns:
str: name of base
"""
......@@ -101,11 +50,11 @@ def getCoordSystemName(base):
return base
def lambdifiedTrafo(baseOld, baseNew):
def lambdifiedTrafo(base_old, base_new):
"""
Args:
baseOld (sympy.CoordSystem)
baseNew (sympy.CoordSystem)
base_old (sympy.CoordSystem)
base_new (sympy.CoordSystem)
Examples:
>>> import numpy as np
......@@ -125,42 +74,73 @@ def lambdifiedTrafo(baseOld, baseNew):
>>> assert new[0, 0] == 5
"""
coords = tuple(baseOld.coord_function(i) for i in range(baseOld.dim))
coords = tuple(base_old.coord_function(i) for i in range(base_old.dim))
f = sympy.lambdify(coords,
baseOld.coord_tuple_transform_to(baseNew,
base_old.coord_tuple_transform_to(base_new,
list(coords)),
modules='numpy')
return f
def transform(array, baseOld, baseNew):
def transform(array, base_old, base_new):
"""
Transform the input array in place
Args:
array (np.ndarray)
base_old (str or sympy.CoordSystem):
base_new (str or sympy.CoordSystem):
Examples:
Cylindrical coordinates
>>> import tfields
>>> cart = np.array([[0, 0, 0],
... [1, 0, 0],
... [1, 1, 0],
... [0, 1, 0],
... [-1, 1, 0],
... [-1, 0, 0],
... [-1, -1, 0],
... [0, -1, 0],
... [1, -1, 0],
... [0, 0, 1]])
>>> cyl = tfields.bases.transform(cart, 'cartesian', 'cylinder')
>>> cyl
Transform cylinder to spherical. No connection is defined so routing via
cartesian
>>> import numpy as np
>>> import tfields
>>> b = np.array([[5, np.arctan(4. / 3), 0]])
>>> newB = tfields.bases.transform(b, 'cylinder', 'spherical')
>>> newB = b.copy()
>>> tfields.bases.transform(b, 'cylinder', 'spherical')
>>> assert newB[0, 0] == 5
>>> assert round(newB[0, 1], 10) == round(b[0, 1], 10)
"""
baseOld = getCoordSystem(baseOld)
baseNew = getCoordSystem(baseNew)
if baseNew not in baseOld.transforms:
for baseTmp in baseNew.transforms:
if baseTmp in baseOld.transforms:
tmpArray = transform(array, baseOld, baseTmp)
return transform(tmpArray, baseTmp, baseNew)
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.")
trafo = tfields.bases.lambdifiedTrafo(baseOld, baseNew)
# very fast trafos in numpy only
shortTrafo = None
try:
shortTrafo = getattr(base_old, 'to_{base_new.name}'.format(**locals()))
except AttributeError:
pass
if shortTrafo:
shortTrafo(array)
return
# trafo via lambdified sympy expressions
trafo = tfields.bases.lambdifiedTrafo(base_old, base_new)
with warnings.catch_warnings():
warnings.filterwarnings('ignore', message="invalid value encountered in double_scalars")
new = np.concatenate([trafo(*coords).T for coords in array])
return new
array[:] = np.concatenate([trafo(*coords).T for coords in array])
if __name__ == '__main__': # pragma: no cover
......
import tfields
import sympy
import numpy as np
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.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)
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)
def cartesian_to_cylinder(array):
x_vals = np.copy(array[:, 0])
y_vals = np.copy(array[:, 1])
problemPhiIndices = 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')
# phi
array[:, 1] = np.sign(y_vals) * np.arccos(x_vals / array[:, 0])
np.seterr(divide='warn', invalid='warn')
array[:, 1][problemPhiIndices] = np.pi
tfields.convert_nan(array, 0.) # for cases like cartesian 0, 0, 1
cylinder.to_cartesian = cylinder_to_cartesian
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) * sympy.sin(theta),
r * sympy.sin(phi) * sympy.sin(theta),
r * sympy.cos(theta)],
inverse=False,
fill_in_gaps=False)
cartesian.connect_to(spherical,
[x, y, z],
[sympy.sqrt(x**2 + y**2 + z**2),
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)
This diff is collapsed.
#!/usr/bin/env
# encoding: utf-8
"""
Author: Daniel Boeckenhoff
Mail: daniel.boeckenhoff@ipp.mpg.de
Collection of additional numpy functions
part of tfields library
"""
import numpy as np
def convert_nan(array, value=0.):
"""
Replace all occuring NaN values by value
"""
nanIndices = np.isnan(array)
array[nanIndices] = value
def view1D(ar):
"""
https://stackoverflow.com/a/44999009/ @Divakar
"""
ar = np.ascontiguousarray(ar)
voidDt = np.dtype((np.void, ar.dtype.itemsize * ar.shape[1]))
return ar.view(voidDt).ravel()
def argsort_unique(idx):
"""
https://stackoverflow.com/a/43411559/ @Divakar
"""
n = idx.size
sidx = np.empty(n, dtype=int)
sidx[idx] = np.arange(n)
return sidx
def duplicates(ar, axis=None):
"""
View1D version of duplicate search
Speed up version after
https://stackoverflow.com/questions/46284660 \
/python-numpy-speed-up-2d-duplicate-search/46294916#46294916
Args:
ar (array_like): array
other args: see np.isclose
Examples:
>>> import tfields
>>> import numpy as np
>>> a = np.array([[1, 0, 0], [1, 0, 0], [2, 3, 4]])
>>> tfields.duplicates(a, axis=0)
array([0, 0, 2])
Returns:
list of int: int is pointing to first occurence of unique value
"""
if axis != 0:
raise NotImplementedError()
sidx = np.lexsort(ar.T)
b = ar[sidx]
groupIndex0 = np.flatnonzero((b[1:] != b[:-1]).any(1)) + 1
groupIndex = np.concatenate(([0], groupIndex0, [b.shape[0]]))
ids = np.repeat(range(len(groupIndex) - 1), np.diff(groupIndex))
sidx_mapped = argsort_unique(sidx)
ids_mapped = ids[sidx_mapped]
grp_minidx = sidx[groupIndex[:-1]]
out = grp_minidx[ids_mapped]
return out
def index(ar, entry, rtol=0, atol=0, equal_nan=False, axis=None):
"""
Examples:
>>> import tfields
>>> a = np.array([[1, 0, 0], [1, 0, 0], [2, 3, 4]])
>>> tfields.index(a, [2, 3, 4], axis=0)
2
>>> a = np.array([[1, 0, 0], [2, 3, 4]])
>>> tfields.index(a, 4)
5
Returns:
list of int: indices of point occuring
"""
if axis is None:
ar = ar.flatten()
elif axis != 0:
raise NotImplementedError()
for i, part in enumerate(ar):
isclose = np.isclose(part, entry, rtol=rtol, atol=atol,
equal_nan=equal_nan)
if axis is not None:
isclose = isclose.all()
if isclose:
return i
if __name__ == '__main__':
import doctest
doctest.testmod()
else:
from . import grid
from .grid import igrid
from . import stats
from .stats import mode, median, mean
import numpy as np
def igrid(*base_vectors, **kwargs):
"""
Args:
*base_vectors (tuple, list or np.array): base vectors spaning the grid
behaviour for different input types:
tuple: will be transformed to slices and given to np.mgrid
list or np.array: arbitrary base vectors
**kwargs
iter_order (list): order in which the iteration will be done.
Frequency rises with position in list. default is [0, 1, 2]
iteration will be done like::
for v0 in base_vectors[iter_order[0]]:
for v1 in base_vectors[iter_order[1]]:
for v2 in base_vectors[iter_order[2]]:
...
Examples:
Initilaize using the mgrid notation
>>> import tfields
>>> tfields.igrid((0, 1, 2j), (3, 4, 2j), (6, 7, 2j))
array([[ 0., 3., 6.],
[ 0., 3., 7.],
[ 0., 4., 6.],
[ 0., 4., 7.],
[ 1., 3., 6.],
[ 1., 3., 7.],
[ 1., 4., 6.],
[ 1., 4., 7.]])
>>> import numpy as np
>>> tfields.igrid([3, 4], np.linspace(0, 1, 2),
... (6, 7, 2), iter_order=[1, 0, 2])
array([[ 3., 0., 6.],
[ 3., 0., 7.],
[ 4., 0., 6.],
[ 4., 0., 7.],
[ 3., 1., 6.],
[ 3., 1., 7.],
[ 4., 1., 6.],
[ 4., 1., 7.]])
>>> tfields.igrid(np.linspace(0, 1, 2), np.linspace(3, 4, 2),
... np.linspace(6, 7, 2), iter_order=[2, 0, 1])
array([[ 0., 3., 6.],
[ 0., 4., 6.],
[ 1., 3., 6.],
[ 1., 4., 6.],
[ 0., 3., 7.],
[ 0., 4., 7.],
[ 1., 3., 7.],
[ 1., 4., 7.]])
"""
iter_order = kwargs.pop('iter_order', np.arange(len(base_vectors)))
# ensure, that the third entry in base_vector of type tuple becomes a complex type
base_vectors = list(base_vectors)
for i, vector in enumerate(base_vectors):
if isinstance(vector, tuple):
if len(vector) == 3:
vector = list(vector)
vector[2] = complex(vector[2])
base_vectors[i] = tuple(vector)
# create the grid
if all([isinstance(val, tuple) for val in base_vectors]):
base_vectors = [slice(*base_vectors[index]) for index in iter_order]
obj = np.mgrid[base_vectors]
obj = obj.reshape(obj.shape[0], -1).T
else:
base_vectors = list(base_vectors)
# transform tuples to arrays with mgrid
for i in range(len(base_vectors)):
if isinstance(base_vectors[i], tuple):
base_vectors[i] = np.mgrid[slice(*base_vectors[i])]
if isinstance(base_vectors[i], list):
base_vectors[i] = np.array(base_vectors[i])
obj = np.empty(shape=(reduce(lambda x, y: x * y,
map(len, base_vectors)),
len(base_vectors)))
def loop_rec(y, n_max, i=0, n=None, *vals):
if n is None:
n = n_max
if n > 0:
for val in y[n_max - n]:
i = loop_rec(y, n_max, i, n - 1, val, *vals)
else:
obj[i] = list(reversed(vals))
i += 1
return i