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

array and grid but still buggy

parent 08eb5691
from .array import igrid
from . import core from . import core
from .array_set_ops import duplicates, index, convert_nan # NOQA
from . import bases from . import bases
from . import array
# __all__ = ['core', 'points3D'] # __all__ = ['core', 'points3D']
from .core import Tensors, TensorFields, TensorMaps from .core import Tensors, TensorFields, TensorMaps
......
...@@ -8,8 +8,8 @@ Collection of additional numpy functions ...@@ -8,8 +8,8 @@ Collection of additional numpy functions
part of tfields library part of tfields library
""" """
import numpy as np import numpy as np
import loggingTools from . import grid
logger = loggingTools.Logger(__name__) from .grid import igrid
def convert_nan(array, value=0.): def convert_nan(array, value=0.):
......
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))
Tensors([[ 0., 3., 6.],
[ 0., 3., 7.],
[ 0., 4., 6.],
[ 0., 4., 7.],
[ 1., 3., 6.],
[ 1., 3., 7.],
[ 1., 4., 6.],
[ 1., 4., 7.]])
>>> tfields.Tensors.grid(np.linspace(3, 4, 2), np.linspace(0, 1, 2),
... np.linspace(6, 7, 2), iter_order=[1, 0, 2])
Tensors([[ 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.Tensors.grid(np.linspace(0, 1, 2), np.linspace(3, 4, 2),
... np.linspace(6, 7, 2), iter_order=[2, 0, 1])
Tensors([[ 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)))
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
loop_rec([base_vectors[i] for i in iter_order], len(base_vectors))
swap_indices = compare_permutations(iter_order, np.arange(len(base_vectors)))
swap_columns(obj, *swap_indices)
return obj
def swap_columns(array, *index_tuples):
"""
Args:
array (list or array)
expects tuples with indices to swap as arguments
Examples:
>>> import tfields
>>> l = np.array([[3,3,3], [2,2,2], [1,1,1], [0, 0, 0]])
>>> tfields.array.grid.swap_columns(l, (1, 2), (0, 3))
>>> l
array([[0, 0, 0], [1, 1, 1], [2, 2, 2], [3, 3, 3]])
"""
# test swap_indices type
for si in index_tuples:
if hasattr(si, '__iter__'):
if len(si) == 2:
continue
raise TypeError("swap_indices must be tuple but is {}"
.format(si))
for i, j in index_tuples:
array[:, [i, j]] = array[:, [j, i]]
def swap_rows(array, *args):
"""
Args:
array (list)
expects tuples with indices to swap as arguments
Examples:
>>> import tfields
>>> l = [[3,3,3], [2,2,2], [1,1,1], [0, 0, 0]]
>>> tfields.array.swap_rows(l, (1, 2), (0, 3))
>>> l
[[0, 0, 0], [1, 1, 1], [2, 2, 2], [3, 3, 3]]
"""
for i, j in args:
array[i], array[j] = array[j], array[i]
def compare_permutations(permut1, permut2):
"""
Return what you need to switch in order to make permut1 become permut2
Examples:
>>> import tfields
>>> a = [1, 2, 0, 4, 3]
>>> b = [0, 1, 2, 3, 4]
>>> si = tfields.array.grid.compare_permutations(a, b)
>>> si
[(0, 2), (1, 2), (3, 4)]
>>> a = np.array(a)
>>> tfields.array.grid.swap_columns(a, *si)
>>> a = list(a)
>>> a
[0, 1, 2, 3, 4]
>>> a == b
True
"""
swap_indices = []
permut1 = list(permut1)
i = 0
while i < len(permut2):
if not permut1[i] == permut2[i]:
j = permut1.index(permut2[i])
swap_rows(permut1, (i, j))
swap_indices.append((i, j))
i += 1
return swap_indices
if __name__ == '__main__':
import doctest
doctest.testmod()
...@@ -11,6 +11,10 @@ import tfields.bases ...@@ -11,6 +11,10 @@ import tfields.bases
import numpy as np import numpy as np
from contextlib import contextmanager from contextlib import contextmanager
from collections import Counter from collections import Counter
import sympy
import scipy as sp
import scipy.spatial # NOQA: F401
import scipy.stats as stats
np.seterr(all='warn', over='raise') np.seterr(all='warn', over='raise')
...@@ -326,6 +330,62 @@ class Tensors(AbstractNdarray): ...@@ -326,6 +330,62 @@ class Tensors(AbstractNdarray):
tensors = np.append(tensors, obj, axis=0) tensors = np.append(tensors, obj, axis=0)
return cls.__new__(cls, tensors, **kwargs) return cls.__new__(cls, tensors, **kwargs)
@classmethod
def grid(cls, *base_vectors, **kwargs):
"""
Args:
baseVector 0 (list/np.array of base coordinates)
baseVector 1 (list/np.array of base coordinates)
baseVector 2 (list/np.array of base coordinates)
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]]:
coords0.append(locals()['v%i' % iter_order[0]])
coords1.append(locals()['v%i' % iter_order[1]])
coords2.append(locals()['v%i' % iter_order[2]])
Examples:
Initilaize using the mgrid notation
>>> npTools.Tensors.grid((0, 1, 2j), (3, 4, 2j), (6, 7, 2j))
Tensors([[ 0., 3., 6.],
[ 0., 3., 7.],
[ 0., 4., 6.],
[ 0., 4., 7.],
[ 1., 3., 6.],
[ 1., 3., 7.],
[ 1., 4., 6.],
[ 1., 4., 7.]])
>>> npTools.Tensors.grid(np.linspace(3, 4, 2), np.linspace(0, 1, 2),
... np.linspace(6, 7, 2), iter_order=[1, 0, 2])
Tensors([[ 3., 0., 6.],
[ 3., 0., 7.],
[ 4., 0., 6.],
[ 4., 0., 7.],
[ 3., 1., 6.],
[ 3., 1., 7.],
[ 4., 1., 6.],
[ 4., 1., 7.]])
>>> npTools.Tensors.grid(np.linspace(0, 1, 2), np.linspace(3, 4, 2),
... np.linspace(6, 7, 2), iter_order=[2, 0, 1])
Tensors([[ 0., 3., 6.],
[ 0., 4., 6.],
[ 1., 3., 6.],
[ 1., 4., 6.],
[ 0., 3., 7.],
[ 0., 4., 7.],
[ 1., 3., 7.],
[ 1., 4., 7.]])
"""
inst = cls.__new__(cls, tfields.array.igrid(*base_vectors, **kwargs))
return inst
@property @property
def rank(self): def rank(self):
""" """
...@@ -432,6 +492,50 @@ class Tensors(AbstractNdarray): ...@@ -432,6 +492,50 @@ class Tensors(AbstractNdarray):
self.transform(baseBefore) self.transform(baseBefore)
def mirror(self, coordinate, condition=None):
"""
Reflect/Mirror the entries meeting <condition> at <coordinate> = 0
Args:
coordinate (int): coordinate index
Examples:
>>> p = tfields.Tensors([[1., 2., 3.], [4., 5., 6.], [1, 2, -6]])
>>> p.mirror(1)
>>> p
Tensors([[ 1., -2., 3.],
[ 4., -5., 6.],
[ 1., -2., -6.]])
multiple coordinates can be mirrored. Eg. a point mirrorion would be
>>> p = tfields.Tensors([[1., 2., 3.], [4., 5., 6.], [1, 2, -6]])
>>> p.mirror([0,2])
>>> p
Tensors([[-1., 2., -3.],
[-4., 5., -6.],
[-1., 2., 6.]])
You can give a condition as mask or as str.
The mirroring will only be applied to the points meeting the condition.
>>> import sympy
>>> x, y, z = sympy.symbols('x y z')
>>> p.mirror([0,2], y > 3)
>>> p
Tensors([[-1., 2., -3.],
[ 4., 5., 6.],
[-1., 2., 6.]])
"""
if condition is None:
condition = np.array([True for i in range(len(self))])
elif isinstance(condition, sympy.Basic):
condition = self.getMask(condition)
if isinstance(coordinate, list) or isinstance(coordinate, tuple):
for c in coordinate:
self.mirror(c, condition)
elif isinstance(coordinate, int):
self[:, coordinate][condition] *= -1
else:
raise TypeError()
def to_segment(self, segment, num_segments, coordinate, def to_segment(self, segment, num_segments, coordinate,
periodicity=2 * np.pi, offset=0, periodicity=2 * np.pi, offset=0,
coordSys=None): coordSys=None):
...@@ -503,6 +607,32 @@ class Tensors(AbstractNdarray): ...@@ -503,6 +607,32 @@ class Tensors(AbstractNdarray):
""" """
return any(self.equal(other, return_bool=False).all(1)) return any(self.equal(other, return_bool=False).all(1))
def indices(self, tensor):
"""
Returns:
list of int: indices of tensor occuring
"""
indices = []
for i, p in enumerate(self):
if all(p == tensor):
indices.append(i)
return indices
def index(self, tensor):
"""
Args:
tensor
Returns:
int: index of tensor occuring
"""
indices = self.indices(tensor)
if not indices:
return None
if len(indices) == 1:
return indices[0]
raise ValueError("Multiple occurences of value {}"
.format(tensor))
def getMoment(self, moment): def getMoment(self, moment):
""" """
Returns: Returns:
...@@ -526,6 +656,185 @@ class Tensors(AbstractNdarray): ...@@ -526,6 +656,185 @@ class Tensors(AbstractNdarray):
else: else:
raise NotImplementedError("Moment %i not implemented." % moment) raise NotImplementedError("Moment %i not implemented." % moment)
def closestPoints(self, other, **kwargs):
"""
Args:
other (Tensors): closest points to what? -> other
**kwargs: forwarded to scipy.spatial.cKDTree.query
Returns:
array shape(len(self)): Indices of other points that are closest to own points
Examples:
>>> m = tfields.Tensors([[1,0,0], [0,1,0], [1,1,0], [0,0,1],
... [1,0,1]])
>>> p = tfields.Tensors([[1.1,1,0], [0,0.1,1], [1,0,1.1]])
>>> p.closestPoints(m)
array([2, 3, 4])
"""
with other.tmp_transform(self.coordSys):
# balanced_tree option gives huge speedup!
kdTree = sp.spatial.cKDTree(other, 1000,
balanced_tree=False)
res = kdTree.query(self, **kwargs)
array = res[1]
return array
def getMask(self, cutExpression=None, coordSys=None):
"""
Args:
cutExpression (sympy logical expression)
coordSys (str): coordSys to evaluate the expression in.
Returns: np.array of dtype bool with lenght of number of points in self.
This array is True, where cutExpression evaluates True.
Examples:
>>> x, y, z = sympy.symbols('x y z')
>>> p = tfields.Points3D([[1., 2., 3.], [4., 5., 6.], [1, 2, -6],
... [-5, -5, -5], [1,0,-1], [0,1,-1]])
>>> p.getMask(x > 0)
array([ True, True, True, False, True, False], dtype=bool)
And combination
>>> p.getMask((x > 0) & (y < 3))
array([ True, False, True, False, True, False], dtype=bool)
Or combination
>>> p.getMask((x > 0) | (y > 3))
array([ True, True, True, False, True, False], dtype=bool)
"""
coords = sympy.symbols('x y z')
with self.tmp_transform(coordSys):
mask = tfields.getMask(self, cutExpression, coords=coords)
return mask
def cut(self, cutExpression, coordSys=None):
"""
Default cut method for Points3D. Works on a copy.
Args:
cutExpression (sympy logical expression): logical expression which will be evaluated.
use symbols x, y and z
coordSys (str): coordSys to evaluate the expression in.
Examples:
>>> x, y, z = sympy.symbols('x y z')
>>> p = tfields.Points3D([[1., 2., 3.], [4., 5., 6.], [1, 2, -6],
... [-5, -5, -5], [1,0,-1], [0,1,-1]])
>>> p.cut(x > 0)
Points3D([[ 1., 2., 3.],
[ 4., 5., 6.],
[ 1., 2., -6.],
[ 1., 0., -1.]])
combinations of cuts
>>> p.cut((x > 0) & (z < 0))
Points3D([[ 1., 2., -6.],
[ 1., 0., -1.]])
Returns:
copy of self with cut applied
"""
if len(self) == 0:
return self.copy()
mask = self.getMask(cutExpression, coordSys=coordSys)
mask.astype(bool)
inst = self[mask].copy()
return inst
def mapToCut(self, cutOrMeshMap, coordSys=None, atIntersection="remove"):
"""
Args:
cutOrMeshMap (cutExpression / Mesh3D with faceScalars)
"""
# in current implementation Points3D use always cut method, never MeshMap
return self.cut(cutOrMeshMap, coordSys=coordSys), None
def to3DArrays(self, arrayShape):
"""
Args:
arrayShape (tuple of lenght 3): how coordinates will be devided.
If you initialize Points3D with a grid of dimensions a,b,c and
want to retrieve the grid-points use (a,b,c)
"""
if not isinstance(arrayShape, tuple):
raise TypeError("ArrayShape is not of type tuple but %s" % type(arrayShape))
if not len(arrayShape) == 3:
raise TypeError("ArrayShape is not of length 3 but %i" % len(arrayShape))
return self.T.reshape((3,) + arrayShape)
def distances(self, other, metric='euclidean', **kwargs):
"""
Examples:
>>> p = tfields.Tensors.grid((0, 2, 3),
... (0, 2, 3),
... (0, 0, 1))
>>> p[4,2] = 1
>>> p.distances(p)[0,0]
0.0
>>> p.distances(p)[5,1]
1.4142135623730951
>>> p.distances([[0,1,2]])[-1]
array([ 3.])
"""
return sp.spatial.distance.cdist(self, other, metric=metric, **kwargs)
def minDistances(self, other=None, metric='euclidean', **kwargs):
"""
Examples:
>>> p = tfields.Tensors.grid((0, 2, 3),
... (0, 2, 3),
... (0, 0, 1))
>>> p[4,2] = 1
>>> dMin = p.minDistances()
>>> dMin
array([ 1. , 1. , 1. , 1. , 1.41421356,
1. , 1. , 1. , 1. ])
>>> dMin2 = p.minDistances(memorySaving=True)
>>> bool((dMin2 == dMin).all())
True
"""
memorySaving = kwargs.pop('memorySaving', False)
if other is not None:
raise NotImplementedError("Should be easy but make shure not to remove 0s")
else:
other = self
try:
if memorySaving:
raise MemoryError()
d = self.distances(other, metric=metric, **kwargs)
return d[d > 0].reshape(d.shape[0], - 1).min(axis=1)
except MemoryError:
minDists = np.empty(self.shape[0])
for i, point in enumerate(other):
d = self.distances([point], metric=metric, **kwargs)
minDists[i] = d[d > 0].reshape(-1).min()
return minDists
def epsilonNeighbourhood(self, epsilon):
"""
Returns:
indices for those sets of points that lie within epsilon around the other
Examples:
Create mesh grid with one extra point that will have 8 neighbours
within epsilon
>>> p = tfields.Tensors.grid((0, 1, 2j),
... (0, 1, 2j),
... (0, 1, 2j))
>>> p = tfields.Tensors([p, [0.5, 0.5, 0.5]])
>>> [len(en) for en in p.epsilonNeighbourhood(0.9)]
[2, 2, 2, 2, 2, 2, 2, 2, 9]
"""
indices = np.arange(self.shape[0])
dists = self.distances(self)
distsInEpsilon = dists <= epsilon
return [indices[die] for die in distsInEpsilon]
......
...@@ -13,10 +13,6 @@ import osTools ...@@ -13,10 +13,6 @@ import osTools
import ioTools import ioTools
import pyTools import pyTools
import mplTools as mpt import mplTools as mpt
import sympy
import scipy as sp
import scipy.spatial # NOQA: F401
import scipy.stats as stats
import warnings import warnings
import loggingTools import loggingTools
logger = loggingTools.Logger(__name__) logger = loggingTools.Logger(__name__)
...@@ -145,9 +141,6 @@ class Points3D(tfields.Tensors): ...@@ -145,9 +141,6 @@ class Points3D(tfields.Tensors):
[-5.000000e+00, -5.000000e+00, -5.000000e+00], [-5.000000e+00, -5.000000e+00, -5.000000e+00],
[ 1.000000e+00, 0.000000e+00, -1.000000e+00], [ 1.000000e+00, 0.000000e+00, -1.000000e+00],
[ 6.123234e-17, 1.000000e+00, -1.000000e+00]]) [ 6.123234e-17, 1.000000e+00, -1.000000e+00]])
>>> p.getPhi()*180/3.1415
Points3D([ 63.43681974, 51.34170594, 63.43681974, -135.00398161,
0. , 90.00265441])
Appending to the array with numpy.append Appending to the array with numpy.append
>>> a = tfields.Points3D([[1,2,3]]) >>> a = tfields.Points3D([[1,2,3]])
...@@ -163,18 +156,6 @@ class Points3D(tfields.Tensors): ...@@ -163,18 +156,6 @@ class Points3D(tfields.Tensors):
kwargs['dim'] = 3 kwargs['dim'] = 3
return super(Points3D, cls).__new__(cls, tensors, **kwargs) return super(Points3D, cls).__new__(cls, tensors, **kwargs)
def contains(self, other):
"""
Inspired by a speed argument @
stackoverflow.com/questions/14766194/testing-whether-a-numpy-array-contains-a-given-row
Examples:
>>> p = tfields.Points3D([[1,2,3], [4,5,6], [6,7,8]])
>>> p.contains([4,5,6])
True
"""
return any(np.equal(self, other).all(1))
def save(self, filePath, extension='npz', **kwargs): def save(self, filePath, extension='npz', **kwargs):
""" """
save a tensors object. save a tensors object.
...@@ -372,250 +353,6 @@ class Points3D(tfields.Tensors): ...@@ -372,250 +353,6 @@ class Points3D(tfields.Tensors):
raise NotImplementedError("Do not understand the format of file {0}".format(filePath)) raise NotImplementedError("Do not understand the format of file {0}".format(filePath))
return fun(filePath, *args, **kwargs) return fun(filePath, *args, **kwargs)
def indices(self, point):
"""
Returns:
list of int: indices of point occuring
"""
indices = []
for i, p in enumerate(self):
if all(p == point):
indices.append(i)
return indices
def index(self, point):
"""
Args: