Skip to content

GitLab

Commit 95b91135 authored by Daniel Boeckenhoff's avatar Daniel Boeckenhoff
Browse files

all doctests run!

parent 164e657b
Expand all Hide whitespace changes
Inline Side-by-side
tfields/__init__.py
View file @ 95b91135
......@@ -17,6 +17,4 @@ from .points3D import Points3D # NOQA
from .mesh3D import Mesh3D # NOQA
from .mesh3D import fields_to_scalars, scalars_to_fields
from .triangles3D import Triangles3D # NOQA
from .scalarField3D import ScalarField3D # NOQA
from .vectorField3D import VectorField3D # NOQA
from .planes3D import Planes3D # NOQA
tfields/bases/__init__.py
View file @ 95b91135
......@@ -7,146 +7,7 @@ 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 .bases import *
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.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):
bse_tpe = type(base)
expctd_tpe = type(sympy.diffgeom.CoordSystem)
raise TypeError("Wrong type of coordSystem base {bse_tpe}. "
"Expected {expctd_tpe}"
.format(**locals()))
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
>>> 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 = b.copy()
>>> 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()
tfields/bases/bases.py 0 → 100644
View file @ 95b91135
#!/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
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):
bse_tpe = type(base)
expctd_tpe = type(sympy.diffgeom.CoordSystem)
raise TypeError("Wrong type of coordSystem base {bse_tpe}. "
"Expected {expctd_tpe}"
.format(**locals()))
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
>>> 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 = b.copy()
>>> 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()
tfields/core.py deleted 100644 → 0
View file @ 164e657b
This diff is collapsed.
tfields/lib/sets.py 0 → 100644
View file @ 95b91135
"""
Algorithms around set operations
"""
import numpy as np
import tfields
class UnionFind(object):
"""
Source:
http://code.activestate.com/recipes/215912-union-find-data-structure/
This algorithm and data structure are primarily used for Kruskal's Minimum
Spanning Tree algorithm for graphs, but other uses have been found.
The Union Find data structure is not a universal set implementation, but can
tell you if two objects are in the same set, in different sets, or you can
combine two sets.
ufset.find(obja) == ufset.find(objb)
ufset.find(obja) != ufset.find(objb)
ufset.union(obja, objb)
"""
def __init__(self):
"""
Create an empty union find data structure.
"""
self.num_weights = {}
self.parent_pointers = {}
self.num_to_objects = {}
self.objects_to_num = {}
self.__repr__ = self.__str__
def insert_objects(self, objects):
"""
Insert a sequence of objects into the structure. All must be Python hashable.
"""
for obj in objects:
self.find(obj)
def find(self, obj):
"""
Find the root of the set that an object 'obj' is in.
If the object was not known, will make it known, and it becomes its own set.
Object must be Python hashable.'''
"""
if obj not in self.objects_to_num:
obj_num = len(self.objects_to_num)
self.num_weights[obj_num] = 1
self.objects_to_num[obj] = obj_num
self.num_to_objects[obj_num] = obj
self.parent_pointers[obj_num] = obj_num
return obj
stk = [self.objects_to_num[obj]]
par = self.parent_pointers[stk[-1]]
while par != stk[-1]:
stk.append(par)
par = self.parent_pointers[par]
for i in stk:
self.parent_pointers[i] = par
return self.num_to_objects[par]
def union(self, object1, object2):
"""
Combine the sets that contain the two objects given.
Both objects must be Python hashable.
If either or both objects are unknown, will make them known, and combine them.
"""
o1p = self.find(object1)
o2p = self.find(object2)
if o1p != o2p:
on1 = self.objects_to_num[o1p]
on2 = self.objects_to_num[o2p]
w1 = self.num_weights[on1]
w2 = self.num_weights[on2]
if w1 < w2:
o1p, o2p, on1, on2, w1, w2 = o2p, o1p, on2, on1, w2, w1
self.num_weights[on1] = w1 + w2
del self.num_weights[on2]
self.parent_pointers[on2] = on1
def __str__(self):
"""
Included for testing purposes only.
All information needed from the union find data structure can be attained
using find.
"""
sets = {}
for i in xrange(len(self.objects_to_num)):
sets[i] = []
for i in self.objects_to_num:
sets[self.objects_to_num[self.find(i)]].append(i)
out = []
for i in sets.itervalues():
if i:
out.append(repr(i))
return ', '.join(out)
def __call__(self, iterator):
"""
Build unions for whole iterators of any size
"""
self.insert_objects(tfields.lib.util.flatten(iterator))
for item in iterator:
for i1, i2 in tfields.lib.util.pairwise(item):
self.union(i1, i2)
def groups(self, iterator):
"""
Return full groups from iterator
"""
groups = {}
keys = []
for item in iterator:
key = self.find(item[0])
if key not in keys:
keys.append(key)
if key not in groups:
groups[key] = []
groups[key].append(item)
return [groups[k] for k in keys]
def group_indices(self, iterator):
"""
Return full groups from iterator
"""
group_indices = {}
keys = []
for i, item in enumerate(iterator):
key = self.find(item[0])
if key not in keys:
keys.append(key)
if key not in group_indices:
group_indices[key] = []
group_indices[key].append(i)
return [group_indices[k] for k in keys]
def disjoint_groups(iterator):
"""
Disjoint groups implementation
Examples:
>>> import tfields
>>> tfields.lib.sets.disjoint_groups([[0, 0, 0, 'A'], [1, 2, 3], [3, 0],
... [4, 4, 4], [5, 4], ['X', 0.42]])
[[[0, 0, 0, 'A'], [1, 2, 3], [3, 0]], [[4, 4, 4], [5, 4]], [['X', 0.42]]]
>>> tfields.lib.sets.disjoint_groups([[0], [1], [2], [3], [0, 1], [1, 2], [3, 0]])
[[[0], [1], [2], [3], [0, 1], [1, 2], [3, 0]]]
Returns:
list: iterator items grouped in disjoint sets
"""
uf = UnionFind()
uf(iterator)
return uf.groups(iterator)
def disjoint_group_indices(iterator):
"""
Examples:
>>> import tfields
>>> tfields.lib.sets.disjoint_group_indices([[0, 0, 0, 'A'], [1, 2, 3],
... [3, 0], [4, 4, 4], [5, 4], ['X', 0.42]])
[[0, 1, 2], [3, 4], [5]]
>>> tfields.lib.sets.disjoint_group_indices([[0], [1], [2], [3], [0, 1], [1, 2], [3, 0]])
[[0, 1, 2, 3, 4, 5, 6]]
Returns:
list: indices of iterator items grouped in disjoint sets
"""
uf = UnionFind()
if isinstance(iterator, np.ndarray):
iterator = iterator.tolist()
uf(iterator)
return uf.group_indices(iterator)
if __name__ == '__main__':
import doctest
doctest.testmod()
tfields/lib/stats.py
View file @ 95b91135
......@@ -18,6 +18,8 @@ def mode(array, axis=0, bins='auto', range=None):
generalisation of the scipy.stats.mode function for floats with binning
Examples:
Forwarding usage:
>>> import tfields
>>> import numpy as np
>>> tfields.lib.stats.mode([[2,2,3], [4,5,3]])
array([[2, 2, 3]])
>>> tfields.lib.stats.mode([[2,2,3], [4,5,3]], axis=1)
......
tfields/lib/symbolics.py 0 → 100644
View file @ 95b91135
"""
sympy helper functions
"""
import sympy
from sympy.logic.boolalg import BooleanFunction
def split_expression(expr):
"""
Return the expression split up in the basic boolean functions.
"""
if isinstance(expr, sympy.Not):
expr = expr.to_nnf()
if isinstance(expr, BooleanFunction):
expressions = expr.args
else:
expressions = [expr]
return expressions
def to_planes(expr):
"""
Resolve BooleanFunctions to retrieve multiple planes
Examples:
>>> import sympy
>>> from tfields.lib.symbolics import to_planes, to_plane
>>> x, y, z = sympy.symbols('x y z')
>>> eq1 = 2*x > 0
>>> eq2 = 2*y + 3*z <= 4
>>> p12 = to_planes(eq1 & eq2)
>>> p1 = to_plane(eq1)
>>> p12[0] == p1
True
>>> p2 = to_plane(eq2)
>>> p12[1] == p2
True
"""
expressions = split_expression(expr)
return [to_plane(e) for e in expressions]
def to_plane(expr):
"""
Tranlate the expression (coordinate form) to normal form and return as Plane
Examples:
Get 3-d plane for linear equations
>>> import sympy
>>> from tfields.lib.symbolics import to_plane
>>> x, y, z = sympy.symbols('x y z')
>>> eq1 = 2*x - 4
>>> p1 = to_plane(eq1)
>>> assert eq1, p1.equation()
# multiple dimensions work
>>> eq2 = x + 2*y + 3*z - 4
>>> p2 = to_plane(eq2)
>>> assert eq2, p2.equation()
The base point is calculated independent of the coords
>>> eq3 = 2*y + 3*z - 4
>>> p3 = to_plane(eq3)
>>> assert eq3, p3.equation()
Inequalities will be treated like equations
>>> ie1 = 2*y + 3*z > 4
>>> p4 = to_plane(ie1)
>>> assert ie1.lhs - ie1.rhs, p4.equation()
Returns: sympy.Plane
"""
x, y, z = sympy.symbols('x y z')
# convert inequalities
if "Than" in str(expr.__class__):
expr = expr.rhs - expr.lhs
poly = expr.as_poly(x, y, z, domain='RR')
if poly is None:
raise ValueError("Expression {expr} can not be described as a polygon."
.format(**locals()))
coords = (x, y, z)
degrees = tuple(poly.degree(c) for c in coords)
if any([d > 1 for d in degrees]):
raise TypeError("Expression {0} can not be represented as plane"
"since the poly dimension is > 1 for one variable!"
.format(expr))
norm = tuple(expr.coeff(coord) for coord in coords)
# now we need any point on the plane
# find the variable to solve for
reduced_expr = expr.copy()
non_zero_coord = None
for coord, degree in zip(coords, degrees):
if degree == 0 or non_zero_coord is not None:
if degree > 0:
reduced_expr = reduced_expr.subs(coord, 0.)
else:
non_zero_coord = coord
non_zero_value = sympy.solve(reduced_expr, non_zero_coord)[0]
point = sympy.Point3D(*[non_zero_value if non_zero_coord == coord else 0.
for coord in coords])
return sympy.Plane(point, normal_vector=norm)
tfields/lib/util.py 0 → 100644
View file @ 95b91135
"""
Various utility functions
"""
import itertools
from six import string_types
def pairwise(iterable):
"""
iterator s -> (s0,s1), (s1,s2), (s2, s3), ...
Source:
https://stackoverflow.com/questions/5434891/iterate-a-list-as-pair-current-next-in-python
Returns:
two iterators, one ahead of the other
"""
a, b = itertools.tee(iterable)
next(b, None)
return zip(a, b)