Commit 08eb5691 authored by Daniel Böckenhoff (Laptop)'s avatar Daniel Böckenhoff (Laptop)
Browse files

added array_Set_ops, bases as file

parent 31e874d3
#!/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
import loggingTools
logger = loggingTools.Logger(__name__)
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()
#!/usr/bin/env
# encoding: utf-8
"""
Author: Daniel Boeckenhoff
Mail: daniel.boeckenhoff@ipp.mpg.de
part of tfields library
Tools for sympy coordinate transformation
"""
import tfields
import numpy as np
import sympy
import sympy.diffgeom
from six import string_types
import warnings
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.get_coord_system)
Return:
sympy.diffgeom.get_coord_system
"""
if isinstance(base, string_types):
base = getattr(tfields.bases, base)
if not isinstance(base, sympy.diffgeom.CoordSystem):
raise TypeError("Wrong type of coordSystem base.")
return base
def get_coord_system_name(base):
"""
Args:
base (str or sympy.diffgeom.get_coord_system)
Returns:
str: name of base
"""
if isinstance(base, sympy.diffgeom.CoordSystem):
base = str(getattr(base, 'name'))
# if not (isinstance(base, string_types) or base is None):
# baseType = type(base)
# raise ValueError("Coordinate system must be string_type."
# " Retrieved value '{base}' of type {baseType}."
# .format(**locals()))
return base
def lambdifiedTrafo(base_old, base_new):
"""
Args:
base_old (sympy.CoordSystem)
base_new (sympy.CoordSystem)
Examples:
>>> import numpy as np
>>> import tfields
Transform cartestian to cylinder or spherical
>>> a = np.array([[3,4,0]])
>>> trafo = tfields.bases.lambdifiedTrafo(tfields.bases.cartesian,
... tfields.bases.cylinder)
>>> new = np.concatenate([trafo(*coords).T for coords in a])
>>> assert new[0, 0] == 5
>>> trafo = tfields.bases.lambdifiedTrafo(tfields.bases.cartesian,
... tfields.bases.spherical)
>>> new = np.concatenate([trafo(*coords).T for coords in a])
>>> assert new[0, 0] == 5
"""
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')
return f
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
>>> print "TEST"
>>> 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
>>> import npTools
>>> pnp = npTools.Points3D(cart)
>>> pnp.transform('cylinder')
>>> pnp
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')
>>> assert newB[0, 0] == 5
>>> assert round(newB[0, 1], 10) == round(b[0, 1], 10)
"""
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
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")
array[:] = np.concatenate([trafo(*coords).T for coords in array])
if __name__ == '__main__': # pragma: no cover
import doctest
doctest.testmod()
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.
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment