#!/usr/bin/env python
import numpy as np
try: from . import vector
except: import vector #Allows it to be run from package dir.
try: from . import simple_cl
except: import simple_cl #Allows it to be run from package dir.
try: from . import mesh
except: import mesh #Allows it to be run from package dir.
import warnings
import os, json
import fractions
from scipy import interpolate
DERIVED_DEBUG = False
if DERIVED_DEBUG:
import traceback
class ShapeError(Exception):
pass
class PathError(Exception):
pass
class ShapeWarning(Warning):
pass
class ToDoError(Exception):
pass
def shape_str(shape):
return '(' + ', '.join(map(lambda d: '-any-' if d < 0 else str(d), shape)) + ')'
def enforce_shape(arr, shape):
arr = np.asarray(arr)
if arr.shape == shape: #Shortcut
return arr
elif arr.shape == (0,):
shape = tuple(map(lambda n: max(n, 0), shape))
return arr.reshape(shape) #If unshaped, just let numpy do it.
elif arr.ndim > len(shape):
while arr.ndim > len(shape):
if arr.shape[-1] == 1:
arr.shape = arr.shape[:-1]
else:
raise ShapeError('input has two many dimensions (%s) to cast into shape %s' % (arr.ndim, shape_str(shape)))
return arr
elif arr.ndim < len(shape):
arr.shape = (1, ) * (len(shape) - arr.ndim) + arr.shape
new_shape = tuple(so if st<=0 else st for so, st in zip(arr.shape, shape))
if arr.shape == new_shape: #Shortcut
return arr
else:
tile = []
for so, st in zip(arr.shape, new_shape):
if so == st: tile.append(1)
elif so == 1: tile.append(st)
else: raise ShapeError('input not castable into array of shape %s' % shape_str(shape))
return np.tile(arr, tile)
def cumsum_exc(s):
cs = np.zeros_like(s)
cs[1:] = np.cumsum(s, axis=0)[:-1]
return cs
CL = simple_cl.CLSession(device='cpu', use_doubles=False)
CODE = os.path.join(os.path.split(__file__)[0], 'cl_src', 'path.cl')
def _check_code(cl=CL):
if 'path_dummy_kernel' not in getattr(cl, 'kernels', {}):
cl.compile_file(CODE)
[docs]def set_opencl_params(use_doubles=None, device=None, context=None):
'''Set the parameters of the OpenCL Session used for computations.
See ``simple_cl`` documentation for details.
*Note:* this function must be called *before* any functions that use it,
or an error will be raised.
Parameters
----------
use_doubles : bool, optional (default: False)
device : string or list, optional (default: 'cpu')
context : pyopencl Context or CLSession object, optional
If specified, causes CLSession to share a context.
'''
if use_doubles is not None: CL.use_doubles = use_doubles
if device is not None: CL.device = device
if context is not None: CL.initialize(context=context)
#From Donnelly, R. J. (2009). Vortex rings in classical and quantum systems.
# Fluid Dynamics Research, 41(5), 051401 1-31.
CORE_MODELS = {
'solid': dict(alpha=7./4, beta=1./4),
'hollow_cv': dict(alpha=2., beta=1./2),
'hollow_cp': dict(alpha=3./2, beta=1./2),
'hollow_st': dict(alpha=1., beta=0.),
'nlse': dict(alpha=1.615, beta=0.615),
'viscous': dict(alpha=2.05, beta=0.558),
}
DEFAULT_CORE_MODEL = 'solid'
FOURTH_ORDER_POLY_COEFF = np.array([
[1., 0., 0., 0., 0.],
[-25./3, +16., -12., +16./3, -1.],
[+70./3, -208/3., +76., -112./3, 22./3],
[-80./3, +96., -128., +224./3, -16.],
[32./3, -128./3, +64., -128./3, +32./3]
])
FOURTH_ORDER_POLY_SUM = np.array([7, 32, 12, 32, 7]) / 90.
[docs]def set_default_core_model(model):
'''Set the default core model used for energy and velocity calculations.
All models adapted from Donnelly, R. J.,
Vortex rings in classical and quantum systems,
Fluid Dynamics Research, 41(5), 051401 1-31 (2009).
Valid options:
* ``"solid"``: (default) solid rotating core, constant volume
* ``"hollow_cv"``: hollow core, constant volume
* ``"hollow_cp"``: hollow core, constant pressure
* ``"hollow_st"``: hollow core with surface tension
* ``"nlse"``: Nonlinear schrodinger equation, i.e. Gross-Pitaevskii equation
* ``"viscous"``: Viscous relaxed core (i.e. Gaussian)
The default is "solid".
'''
global DEFAULT_CORE_MODEL
if model not in CORE_MODELS:
raise ValueError('invalid core model: "%s"\nvalid options: %s' % (model, CORE_MODELS.keys()))
DEFAULT_CORE_MODEL = model
[docs]class Path(object):
'''Object which describes a 3D path, optionally with attached attributes.
Automatically handles the calculation of geometric parameters like tangents,
etc. as special attributes:
* ``tangent``: the tangent vector, calculated as the normalized path derivative
* ``seg_length``: the segment length
* ``L``: the total length, calculated as ``sum(seg_length)``
* ``s``: path length coordinate, defined from the first point in the path
* ``fs_normal``, ``fs_binormal``: the Frenet-Serret (bi-)normal
* ``torsion``, ``curvature``: the torsion and curvature, from the Frenet-Serret basis
* ``s_polygon``, ``L_polygon``: like ``s`` and ``L``, but for a polygonal path
These parameters will appear as attached as attributes once they are created.
If the path is modified in place these derived attributes should be deleted
with :meth:`clear_derived`. (This happens automatically unless the path
array is modified in place.)
Also allows the attachment of arbitrary parameters to the path as arrays
with shape (N, -any-), where N is the number of points in the path.
These can either be defined at creation, or attached as attributes later.
'''
#Derivation procedes in stages; this way a attribute only depends on those
# at lower levels.
derived_attributes = {
#Base attributes
'delta': '_derive_stage1',
'seg_polygon': '_derive_stage1',
'L_polygon': '_derive_stage1',
#Depend on stage 1
's_polygon': '_derive_s_polygon',
#Create weights for taking derivatives
'ds_weights': '_derive_stage2',
#Depends on stage 2
'tangent': '_derive_tangent',
'fs_normal': '_derive_fs',
'curvature': '_derive_fs',
'fs_binormal': '_derive_fs_binormal',
'torsion': '_derive_torsion',
#Create points for a Cubic Bezier-resampler
'cb_P0': '_derive_stage3',
'cb_P1': '_derive_stage3',
'cb_P2': '_derive_stage3',
'cb_P3': '_derive_stage3',
'cb_ds0': '_derive_stage3',
'cb_ds1': '_derive_stage3',
'cb_ds2': '_derive_stage3',
'cb_ds3': '_derive_stage3',
'cb_ds4': '_derive_stage3',
'seg_length': '_derive_stage3',
#Depnds on stage 3
's': '_derive_s',
'L': '_derive_L',
}
special_attributes = {
'info': dict,
'closed': bool,
'derivative_points': int,
}
def __init__(self, path, closed=True, info={}, derivative_points=5, **kwargs):
'''Create a new path. In addition to the parameters above, aribitrary
keyword arguments may be specified which correspond to attached parameters.
Parameters
----------
path : [N, 3] array
closed : bool (default: True)
info : dictionary
Arbitrary info about the path. Will be saved with the path.
derivative points : int (3 or 5, default 5)
The number of points used for deriviative estimation
**kwargs : [N, ...] arrays
Arbitrary attributes may be attached, provided they are castable to
arrays where the first dimension is the same as the number of
points.
'''
self.path = path
self.info = info
self.closed = closed
self.derivative_points = derivative_points
#self.refine_s = refine_s
object.__setattr__(self, 'attached', {})
#object.__setattr__(self, 'derived', {})
for k, v in kwargs.iteritems():
if k in ['normals', 'binormals', 'tangents'] and k[:-1] not in kwargs:
self.__setattr__(k[:-1], v)
else:
self.__setattr__(k, v)
[docs] def clear_derived(self):
'''Clear all derived attributes.
Should be used if the path is modified in place by external commands.'''
#object.__setattr__(self, 'derived', {})
for k in self.derived_attributes.keys():
#if hasattr(self, k): delattr(self, k) #This is bad! Called hasattr causes it to be calculated!
if k in self.__dict__: delattr(self, k)
def __len__(self):
return len(self.path)
def __getattr__(self, k):
#This only gets called if a derived attribute is requested that has
# not already been computed.
if DERIVED_DEBUG:
print '---- Asked for derivation of attribute "%s" ----' % k
traceback.print_stack()
if k in self.derived_attributes:
#if k != 'delta': raise AttributeError("Tried to get '%s'" % k)
getattr(self, self.derived_attributes[k])() #Get function and call it
if k in self.__dict__:
return getattr(self, k)
else:
raise AttributeError("property '%s' should be derived, but derivation function did not define it\n(this should never happen; there is a bug in the source!)" % k)
else:
raise AttributeError("path does not have property '%s'" % k)
#if hasattr(self, 'attached') and k in self.attached:
# return self.attached[k]
#if k in self.derived_attributes:
# getattr(self, self.derived_attributes[k])()
#
#if k in self.derived: #Will catch just-defined properties above
# return self.derived[k]
#
#elif k == 'dtype':
# return self.path.dtype
#
#else:
# raise AttributeError("path does not have property '%s'" % k)
def __setattr__(self, k, v):
if k in self.special_attributes:
object.__setattr__(self, k, self.special_attributes[k](v))
elif k in self.derived_attributes:
raise AttributeError('%s is a derived attribute, and can not be set dircetly') % k
elif k == 'path': #If there are attached objects, enforce that the path shape stays the same
NP = len(self) if ('attached' in self.__dict__ and self.attached.keys()) else -1
object.__setattr__(self, k, enforce_shape(v, (NP, 3)))
object.__setattr__(self, 'dtype', self.path.dtype)
self.clear_derived()
else:
v = np.asarray(v)
if v.ndim == 1 and (len(v) > 10 or len(v) == len(self)) and len(v) != 0:
warnings.warn('property "%s" attached to a path as a 1D array.\nThis is interpreted as a constant vector (of length %d) on each point!\nTo attach a scalar value to the path, reshape the input array to (-1, 1).' % (k, len(v)), ShapeWarning, stacklevel=2)
self.attached[k] = enforce_shape(v, (len(self), -1))
object.__setattr__(self, k, self.attached[k])
def __delattr__(self, k):
object.__delattr__(self, k)
if k in self.attached: del self.attached[k]
[docs] def midpoint(self):
'''Return the segment midpoints.
If ``closed==False``, will return one less than the number of points.
'''
if self.closed:
return 0.5*(vector.plus(self.path) + self.path)
else:
return 0.5*(self.path[1:] + self.path[:-1])
[docs] def copy(self):
'''Return a copy of the path in which all attached arrays are copied.'''
extra = {}
for k, v in self.attached.iteritems():
extra[k] = v.copy()
for k in self.special_attributes.keys():
if hasattr(self, k): extra[k] = getattr(self, k)
return Path(self.path.copy(), **extra)
[docs] def rot(self, rv):
'''Rotate the path in place.
Parameters
----------
rv : vector ([3] array)
Rotation vector, will rotate counter-clockwise around axis by an
amount given be the length of the vector (in radians).
'''
self.apply_basis(vector.rot(rv))
[docs] def scale(self, scale):
'''Scale the path in place.'''
self.path = self.path * scale # using *= would not clear derived attributes!
[docs] def translate(self, offset):
'''Translate the path in place.'''
self.path += offset #Derived attributes should be unaffected.
#--------------------------------------------------------------------------
#---- Derivation, ordered by dependancy
#--------------------------------------------------------------------------
def _derive_stage1(self):
if self.closed:
object.__setattr__(self, 'delta', vector.plus(self.path) - self.path)
else:
object.__setattr__(self, 'delta', self.path[1:] - self.path[:-1])
object.__setattr__(self, 'seg_polygon', vector.mag1(self.delta))
object.__setattr__(self, 'L_polygon', self.seg_polygon.sum())
def _derive_s_polygon(self):
if self.closed:
object.__setattr__(self, 's_polygon', cumsum_exc(self.seg_polygon))
else:
object.__setattr__(self, 's_polygon', np.zeros((len(self), 1), self.dtype))
self.s_polygon[1:] = np.cumsum(self.seg_polygon, axis=0)
def _Ds(self, offset=1):
#Used internally -- note that this routine does not check for closed
# paths. This means the end points should be handled with care.
if offset > 0:
return (np.roll(self.s_polygon, -offset, axis=0) - self.s_polygon) % self.L_polygon
else:
return -((self.s_polygon - np.roll(self.s_polygon, -offset, axis=0)) % self.L_polygon)
def _derive_stage2(self):
if self.derivative_points == 3:
m = self._Ds(-1)
p = self._Ds(+1)
weights = np.zeros((3, len(self), 1), self.dtype)
if self.closed:
weights[0] = -p / (m * (m-p))
weights[1] = -1/m - 1/p
weights[2] = m / (p * (m-p))
else:
weights[0, 1:-1] = -p[1:-1] / (m[1:-1] * (m[1:-1]-p[1:-1]))
weights[0, -1] = +1/m[-1]
weights[1, 0 ] = -1/p[0]
weights[1, 1:-1] = -1/m[1:-1] - 1/p[1:-1]
weights[1, -1] = -1/m[-1]
weights[2, 0 ] = +1/p[0]
weights[2, 1:-1] = m[1:-1] / (p[1:-1] * (m[1:-1]-p[1:-1]))
elif self.derivative_points == 5:
ds = [self._Ds(offset) if offset else 1 for offset in [-2, -1, 0, 1, 2]]
#The interpolation is derived from Lagrange interpolating polynomials
#Here define a function which sums over values, optionally omitting
# some.
def ld(i, start=0, end=5, sl=slice(None, None)):
num = 0 if i == 2 else 1.
den = 1.
for j in range(start, end):
if j == 2: continue
if i == 2:
num -= 1./ds[j][sl]
else:
if i == j: den *= ds[i][sl]
else:
num *= ds[j][sl]
den *= ds[j][sl] - ds[i][sl]
return num/den
weights = np.zeros((5, len(self), 1), self.dtype)
if self.closed:
for i in range(5): weights[i] = ld(i)
else:
sl = slice(2, -2) #Central points
for i in range(0, 5): weights[i][sl] = ld(i, 0, 5, sl=sl)
sl = 0 #First point
for i in range(2, 5): weights[i][sl] = ld(i, 2, 5, sl=sl)
sl = 1 #Second point
for i in range(1, 5): weights[i][sl] = ld(i, 1, 5, sl=sl)
sl = -2 #Second to last point
for i in range(0, 4): weights[i][sl] = ld(i, 0, 4, sl=sl)
sl = -1 #Last point
for i in range(0, 3): weights[i][sl] = ld(i, 0, 3, sl=sl)
object.__setattr__(self, 'ds_weights', weights)
def _derive_tangent(self, method=''):
object.__setattr__(self, 'tangent', vector.norm(self.partial_s('path')))
#return self.derived['tangent']
def _derive_fs(self):
#Access as [path].fs_normal
dT = self.partial_s(self.tangent)
dT -= self.tangent * vector.dot1(self.tangent, dT) #Enforce orthogonality
object.__setattr__(self, 'curvature', vector.mag1(dT))
object.__setattr__(self, 'fs_normal', dT / self.curvature)
def _derive_fs_binormal(self):
#Access as [path].fs_binormal
object.__setattr__(self, 'fs_binormal', vector.cross(self.tangent, self.fs_normal))
def _derive_torsion(self):
#Access as [path].torsion
#minus sign is customary (http://en.wikipedia.org/wiki/Torsion_of_a_curve)
object.__setattr__(self, 'torsion', -vector.dot1(self.fs_normal, self.partial_s(self.fs_binormal)))
def _derive_stage3(self):
object.__setattr__(self, 'cb_P0', self.path if self.closed else self.path[:-1])
object.__setattr__(self, 'cb_P3', np.roll(self.path, -1, axis=0) if self.closed else self.path[1:])
t0 = self.tangent if self.closed else self.tangent[:-1]
t3 = np.roll(self.tangent, -1, axis=0) if self.closed else self.tangent[1:]
phi = np.arccos(np.clip(vector.dot(t0, t3), -1, 1))
#Let's assume this section is a circle. In this case, we can estimate
# the length of the arc segment as (line length) * (1 + (1-t0.t3)/12)
#The distance from P0 to P1 and P2 to P3 will be 1/3 this length.
s = (1./3) * vector.mag1(self.cb_P3 - self.cb_P0) * (1 + (1 - vector.dot1(t0, t3))/12. )
v01 = (t0) * s
v23 = (t3) * s
object.__setattr__(self, 'cb_P1', self.cb_P0 + v01)
object.__setattr__(self, 'cb_P2', self.cb_P3 - v23)
v12 = self.cb_P2 - self.cb_P1
#We'll approximate s(t) as a 4th order polynomial for each segment.
#We find interpolating points by computing the magnitude of the
# derivatives at equally spaced points along the segment.
#dx(t) = (1-t)^3 v01 + 6 (1-t) t v12 + 3 t^2 v23
ds = np.array([
3*s[:, 0],
vector.mag(27./16*v01 + 9./8*v12 + 3./16*v23),
vector.mag(3./4* v01 + 3./2*v12 + 3./4*v23),
vector.mag(3./16* v01 + 9./8*v12 + 27./16*v23),
3*s[:, 0]
]).T
object.__setattr__(self, 'ds_coeff', (FOURTH_ORDER_POLY_COEFF * ds[:, np.newaxis, :]).sum(-1))
object.__setattr__(self, 'seg_length', (FOURTH_ORDER_POLY_SUM * ds[:, np.newaxis, :]).sum(-1))
def _derive_L(self):
object.__setattr__(self, 'L', self.seg_length.sum())
def _derive_s(self):
if self.closed:
object.__setattr__(self, 's', cumsum_exc(self.seg_length))
else:
object.__setattr__(self, 's', np.zeros((len(self), 1), self.dtype))
self.s[1:] = np.cumsum(self.seg_length, axis=0)
[docs] def partial_s(self, a):
'''Take a partial derivative with respect to path length of any quantity on the tangle.
Parameters
----------
a : string or [N, M] array
If a is a string, an attribute of the path is used. Otherwise the
first dimension of a should match the number of points.
Returns
-------
partial_s : [N, M] array
The partial derivative with respect to path length. A 3 or 5 point
stencil is used (depending on self.deriviative_points).
Fully accounts for non-equal segment lengths and end points for
open paths.
'''
if isinstance(a, basestring):
a = getattr(self, a)
a = enforce_shape(a, (len(self), -1))
d = np.zeros_like(a)
for i, w in enumerate(self.ds_weights):
d += w * np.roll(a, self.derivative_points//2-i, axis=0)
return d
[docs] def mean(self, a):
'''Compute the mean of some variable over the path, accounting for non-uniform segment length.
Parameters
----------
a : string or [N, M] array
If a is a string, an attribute of the path is used. Otherwise the
first dimension of a should match the number of points.
Returns
-------
mean : [M] array
'''
if isinstance(a, basestring):
a = getattr(self, a)
a = enforce_shape(a, (len(self), -1))
dL = self.seg_length_per_point()
return (a * dL).sum(0) / dL.sum()
[docs] def seg_length_per_point(self):
'''Compute the segment length associated with each point.
Returns
-------
seg_L_per_point : [N] array
For closed paths, half of the length of each adjacent segment is
assigned to each point. For open paths, the end points only get
a single segment piece.
'''
seg_L = 0.5 * (self.seg_length + np.roll(self.seg_length, 1, axis=0))
if not self.closed:
seg_L[0] -= 0.5 * self.seg_length[-1]
seg_L[-1] -= 0.5 * self.seg_length[0]
return seg_L
#def integrate_s(self, a):
# '''Integrate a quantity on the tangle with respect to path-length.
#
# Uses a trapezoidal rule with the segment lengths.
#
# *TO DO:* a better integrator method (the unequal spacing makes it a tad
# tricky).
#
# Parameters
# ----------
# a : string or [N, M] array
# If a is a string, an attribute of the path is used. Otherwise the
# first dimension of a should match the number of points.
#
# Returns
# -------
# integral_s : [N, M] array
# The integral with respect to path length.
# '''
# if isinstance(a, basestring):
# a = getattr(self, a)
#
# a = enforce_shape(a, (len(self), -1))
#
[docs] def to_dict(self, omit=('derivative_points', 'refine_s')):
'''Convert all parameters to a dictionary, usually for saving.
Parameters
----------
omit : list or tuple
Parameters not to include in dictionary.
Returns
-------
d : dict
'''
keys = set(['path', 'closed', 'info', 'derivative_points', 'refine_s'])
keys.update(self.attached.keys())
d = {}
for k in keys:
if k in omit or k in self.derived_attributes: continue
elif hasattr(self, k): d[k] = getattr(self, k)
return d
[docs] def save(self, fn):
'''Save the path as a file. The filename is not enforced, but ``.path``
is recommended.
The resulting file is a JSON version of ``self.to_dict``.
'''
d = self.to_dict()
d['__file_type__'] = 'JSON 3D PATH'
save_json(d, fn)
[docs] def keep(self, indices):
'''Keep only the specified point indices in the path.
Parameters
----------
indices : int array
'''
i = np.array(indices)
attached = getattr(self, 'attached', {})
object.__setattr__(self, 'attached', {}) #Otherwise setting the path will raise an error!
#print self.attached.keys()
self.path = self.path[i]
for k, v in attached.iteritems():
setattr(self, k, v[i])
[docs] def decimate(self, ratio):
'''Decimate the path points, preserving end points if not closed.
Also decimates all the attached arrays.
Parameters
----------
ratio : int
'''
ratio = int(ratio)
i = np.arange(0, len(self.path), ratio)
end = len(self.path) - 1
if not self.closed and i[-1] != end:
i = np.concatenate([i, [end]])
attached = getattr(self, 'attached', {})
object.__setattr__(self, 'attached', {}) #Otherwise setting the path will raise an error!
#print self.attached.keys()
self.path = self.path[i]
for k, v in attached.iteritems():
setattr(self, k, v[i])
[docs] def eliminate_close_points(self, min_dist):
'''Eliminate points which are too close.
Parameters
----------
min_dist : float
Distance between points below which to clip.
'''
if min_dist <= 0: return
#ip = [0]
#for i in range(1, len(self.path)):
# if vector.mag(self.path[i] - self.path[ip[-1]]) > min_dist:
# ip.append(i)
#
#if self.closed:
# if vector.mag(self.path[0] - self.path[ip[-1]]) < min_dist:
# ip.pop(-1)
#elif ip[-1] != len(self.path) - 1:
# ip.append(len(self.path) - 1)
#
#ip1 = np.array(ip)
good = self.seg_length > min_dist
#Multiple points in a row may have a short distance -- fix this by
# checking for it explicitly, putting points back in as needed.
il = -2
accum = 0
for i in np.where(~good)[0]:
if i == il + 1:
if accum > min_dist:
good[i] = True
accum = 0
else:
accum = 0
accum += self.seg_length[i]
il = i
ip = np.where(good)[0]
#print ip1
#print ip2
#
#test1 = np.zeros(len(self.path), dtype=bool)
#test1[ip1] = True
#
#test2 = np.zeros(len(self.path), dtype=bool)
#test2[ip2] = True
#print self.seg_length[np.where(test1 != test2)]
#ip = ip1
attached = getattr(self, 'attached', {})
object.__setattr__(self, 'attached', {}) #Otherwise setting the path will raise an error!
self.path = self.path[ip]
for k, v in attached.iteritems():
setattr(self, k, v[ip])
[docs] def reduce_by_curvature(self, points_per_circle=100):
r'''Reduce a path by eliminating points in nearly-straight sections.
Works by computing the angle of each corner, and only keeps points when
the total accumulated angle exceeds a threshold given by the
specified number of points per circle.
Parameters
----------
points_per_circle : int, optional (default: 100)
'''
angle_min = 2*np.pi / points_per_circle
n = vector.norm(self.delta)
angle = np.arccos(np.clip(vector.dot(n, np.roll(n, 1, axis=0)), -1, 1))
#If not closed, we MUST keep the end points!
if not self.closed:
angle[0] = angle_min
angle[-1] = angle_min
keep = []
a_s = 0
for i, a in enumerate(angle):
a_s += a
if a_s >= angle_min:
keep.append(i)
a_s = 0
keep = np.array(keep)
attached = getattr(self, 'attached', {})
object.__setattr__(self, 'attached', {}) #Otherwise setting the path will raise an error!
self.path = self.path[keep]
for k, v in attached.iteritems():
setattr(self, k, v[keep])
[docs] def moment(self):
r'''Calculates the path magnetic moment: :math:`\vec{m} = \frac{1}{2} \oint \vec{r} \times d\vec{r}`.
Uses a discrete approximation: :math:`\vec{m} \approx \frac{1}{2} \sum \left(\vec{r}_{i+1} - \vec{r}_i\right) \times \frac{1}{2}\left(\vec{r}_{i+1} + \vec{r}_i\right)`
'''
if not self.closed:
raise PathError('moment is only defined for closed paths')
return 0.5 * np.cross(self.midpoint(), self.delta).sum(0)
[docs] def orient(self, X, Y=None):
'''Orient tangle so specified vector is in the x-direction.
Parameters
----------
X : vector
Y : vector, optional
These vectors will appear as the new x and y axis, respectively.
If only one is specified, the y axis is chosen semi-arbitrarily.
Note: if x and y are not orthogonal, y will be adjusted so that a
right-handed orthonormal basis is formed (in other words the result
is always a pure rotation, with no mirroring, etc.).
'''
B = vector.normalize_basis([X, Y] if Y is not None else [X])
self.apply_basis(B)
[docs] def apply_basis(self, B):
'''Apply an orthonormal basis transformation.
Parameters
----------
B : basis array (3, 3)
Assumed to be orthonormal, but not checked.
'''
self.path = vector.apply_basis(self.path, B)
#for k, v in self.attached.iteritems():
# if 'normal' in k or k == 'tangent':
# self.attached[k] = vector.apply_basis(v, B)
#
#for k, v in self.derived.iteritems():
# if 'normal' in k or k == 'tangent':
# self.derived[k] = vector.apply_basis(v, B)
[docs] def center(self):
r'''Return the center of the path.
Computed as: :math:`\left<\vec{r}\right> = \frac{\sum \frac{L_i}{2}\left(\vec{r}_{i+1} + \vec{r}_i\right)}{\sum L_i}`
'''
return (self.seg_length * self.midpoint()).sum(0) * (1./self.L)
[docs] def rms_radius(self):
r'''Return the root mean squared radius of the path.
Computed as: :math:`\bar{r} = \sqrt{\frac{\sum L_i \left|\frac{\vec{r}_{i+1} + \vec{r}_i}{2} - \left<\vec{r}\right>\right|^2}{\sum L_i}}`
'''
return np.sqrt((self.seg_length * vector.mag1(self.midpoint() - self.center())**2).sum() * (1./self.L))
[docs] def gen_normal(self, normal):
r'''Fix a normal vector to be unit length and perpendicular to the tangent.
'''
return vector.norm(normal - vector.dot1(self.tangent, normal) * self.tangent)
[docs] def twist(self, normal=None, per_point=False, normalize=True):
r'''Measure the twist of a normal vector about a path.
Computed as :math:`Tw = \frac{1}{2\pi}\oint d\vec{r} \cdot \left(\vec{n} \times \partial_s\vec{n}\right)`
If the normal is the Frenet-Serret normal, the twist becomes the
(integrated) torsion; the sign convention is such that right handed
helices have positive twist.
Parameters
----------
normal : vector array, optional (default: self.normal)
per_point : bool, optional (default: False)
If True, returns the twist rate per point instead of the total.
normalize : bool, optional (default: True)
If True, the normal is passed through :meth:`path.gen_normal` first.
Returns
-------
twist : float or array
The total or per point twist density.
'''
if normal is None: normal = self.normal
else: normal = enforce_shape(normal, (len(self), 3))
if normalize: normal = self.gen_normal(normal)
tau = vector.dot(vector.cross(normal, self.partial_s(normal)), self.tangent) / (2*np.pi)
if per_point: return tau
else: return (tau * self.seg_length_per_point().flat).sum()
[docs] def biot_savart_on_path(self, a, gamma=1.0, core_model=None, alpha=None):
'''Compute the biot-savart flow field on the path.
Uses an exact expression for a polygonal path; neighboring segments are
not included in the expression, and are replaced with a computation for
a circular arc going through three points.
Parameters
----------
a : float
The "core size" used in the cutoff calculation.
gamma : float, optional (default: 1)
The core circulation; multiplies the overall flow.
core_model : str, optional
The core model to use; affects the way the short range cutoff is
handlded. The default is specified by a global variable; see
:func:`set_default_core_model` for more information.
alpha : float, optional
The cutoff correction parameter. Usually determined from the
core_model, but if specified directly overrides this setting.
'''
if not self.closed: raise ValueError('biot-savart fields only valid for closed paths')
raise ValueError('not implemented yet!')
[docs] def colorize(self, c, cmap='jet', clim=None):
'''Color a Path using matplotlib specifications.
Parameters
----------
c : string or array
* If type is string, it is interpreted like a normal matplotlib color
(e.g., ``"r"`` for red, or ``"0.5"`` for 50% gray, or ``"#FF8888"`` for pink).
* If type is array, it should have the same length as the number of points,
and it will be converted to a colormap.
cmap : string or matplotlib colormap (default: 'jet')
The color map to use. Only relevant if ``c`` is an array.
clim : tuple (default: None)
The color limits. If None, the max/min of c are used.
'''
if isinstance(c, basestring):
import matplotlib.colors
self.color = np.tile(matplotlib.colors.colorConverter.to_rgb(c), (len(self.path), 1))
else:
import matplotlib.cm
c = enforce_shape(c, (len(self.path),))
if clim is None: clim = (c.min(), c.max())
self.color = matplotlib.cm.get_cmap(cmap)((c - clim[0]) / (clim[1] - clim[0]))
[docs] def color_by_proximity(self, p, c=None, w=None, rgb_only=True, power=-6, period=None):
'''Color a path by proximity to another path, tangle, or set of colored points.
Parameters
----------
p : Path, Tangle or [N, 3] array
c : [N, 3 or 4] array
The color of each point; required if p is an array, and ignored if
it is a Path or Tangle (in this case the ``color`` attribute is
used).
w : [N] array, optional
Weight array.
rgb_only : bool, optional (default: True)
If True, the alpha component is dropped
power : int, optional (default: -6)
The power of the radial weighting function.
Recommended values are <= -4.
period : None, float or (3,) array_like (optional)
If not None, the period to use for wrapping the proximity
calculation (data assumed to be go from 0-period on each axis).
If one number is specified, the period is assumed to be the same
on all three axes.
'''
if hasattr(p, '_color_by_proximity_points'):
p, c, w = p._color_by_proximity_points()
p = enforce_shape(p, (-1, 3))
if c is None: raise ValueError('if specified as points, a color array must also be specified')
c = enforce_shape(c, (len(p), -1))
if c.shape[-1] == 3: c = np.concatenate([c, np.zeros((len(c), 1))], axis=-1)
elif c.shape[-1] != 4: raise ValueError('last dimension of color array must be 3 or 4')
if w is None: w = np.ones(len(c))
else: w = enforce_shape(w, (len(c),) )
_check_code()
pw = CL.to_device(np.concatenate([p, w.reshape(-1, 1)], axis=-1).astype(CL.real))
c = CL.to_device(c.astype(CL.real))
self._color_by_proximity_inner(pw, c, rgb_only, power)
def _color_by_proximity_points(self):
#Internal use
return self.path, self.color, self.seg_length_per_point()
def _color_by_proximity_inner(self, pw, c, rgb_only, power, period=None):
#Internal use only.
x = CL.to_device(self.path.astype(CL.real))
color = CL.empty((len(self.path), 4))
#__kernel void path_color_by_proximity(
# global REAL* X, int nx,
# global REAL4* pw, global REAL4* c, int nc,
# global REAL4* color, int power
# )
if period is None:
CL.path_color_by_proximity(
x, len(self.path),
pw, c, len(pw),
color, power)
else:
p = np.ones(3, dtype=CL.real)
p *= period
px, py, pz = p
CL.path_color_by_proximity_periodic(
x, len(self.path),
pw, c, len(pw),
color, power,
px, py, pz)
color = color.get()
if rgb_only: self.color = color[..., :3]
else: self.color = color
[docs] def cubic_bezier_resample(self, s, attached_interp='linear'):
'''Resample the path using a cubic bezier.
The resulting path will go through all the original points with a
tangent given by the ``tangent`` attribute.
Any properties on the path will be resample as well.
Parameters
----------
s : array
The target s-coordinates (pathlength) of the new points.
attached_interp : str or int, optional (default: ``"linear"``)
The interpolation order for attached attributes.
Passed directly to ``scipy.interpolate.interp1d``.
(Options: ``"linear"``, ``"nearest"``, ``"zero"``, ``"slinear"``,
``"quadratic"``, ``"cubic"``, or an integer specifying the
spline order)
Returns
-------
path : Path
The interpolated path, with all attached attributes resampled.
'''
N = len(self.path)
if self.closed:
t_from_s = interpolate.interp1d(np.concatenate([self.s[:, 0], [self.L]]), np.arange(N+1), kind='linear', copy=False)
t = t_from_s(s % self.L)
i = np.floor(t).astype('i')
else:
t_from_s = interpolate.interp1d(self.s[:, 0], np.arange(N), kind='linear', copy=False)
t = t_from_s(s) #This may trigger a bounds error: we want this!
i = np.floor(t).astype('i')
#Last point gets shifted down one. This point will then have a t
# value of 1, which is ok.
i[np.where(i == len(self.path)-1)] -= 1
t = (t - i)
#We have chosen t assuming the points are equally spaced along the
# bezier resample, which is generally not true. Here we use Newton's
# method to refine the point. Since we started close, 3 cycles should
# be plenty.
#s(t) for each bezier segment is not possible to calculate exactly, so
# we use a 4th order polynomial approximation for which the weights
# are already calculated (in _derive_stage3)
s_target = t * self.seg_length[i, 0]
int_coeff = self.ds_coeff[i, :] * (1./(1+np.arange(5)))
for n in range(3):
s_actual = (int_coeff * t[..., np.newaxis] ** (1+np.arange(5))).sum(-1)
ds_dt = (self.ds_coeff[i, :] * t[..., np.newaxis]**np.arange(5)).sum(-1)
t -= (s_actual - s_target) / ds_dt
#Create new path with our optimized t's
t = t[:, np.newaxis]
ot = 1-t
path = (ot**3 * self.cb_P0[i] +
3*ot**2*t * self.cb_P1[i] +
3*ot*t**2 * self.cb_P2[i] +
t**3 * self.cb_P3[i] )
#Interpolate attached parameters
extra = {}
ep = 1 if attached_interp=='linear' else 4
if self.attached.keys():
if self.closed:
interp_s = np.concatenate([self.s[-ep:, 0]-self.L, self.s[:, 0], self.s[:ep, 0]+self.L])
else:
interp_s = self.s[:, 0]
for k, v in self.attached.iteritems():
if self.closed:
v = np.concatenate([v[-ep:], v, v[:ep]])
extra[k] = interpolate.interp1d(interp_s, v, axis=0, kind=attached_interp, copy=False)(s)
return Path(path, closed=self.closed, info=self.info, **extra)
[docs] def resample_with_spacing(self, dl, attached_interp='linear'):
'''Resample path evenly, spacing points by specified distance.
Parameters
----------
dl : float
New point spacing target. Actual spacing will slightly different,
as the points will be equally spaced.
attached_interp : int or string (default: ``"linear"``)
Type of interpolation for attached attributes.
See :meth:`cubic_bezier_resample` for details.
Returns
-------
path : Path
The interpolated path, with all attached attributes resampled.
'''
s = np.linspace(0, self.L, int(round(self.L / dl)))
if self.closed: s = s[:-1]
return self.cubic_bezier_resample(s, attached_interp=attached_interp)
[docs] def parallel_transport_framing(self, N0=(1, 0, 0), twist_to_close=False):
'''Generate a parallel transport framing, starting from the first point.
Parameters
----------
N0 : [3] array, optional (default: (1, 0, 0))
The normal vector to begin the framing.
twist_to_close : bool, optional (default: False)
If True, the framing is uniformly twisted to close on itself.
Ignored if the path is not closed.
Returns
-------
normal : [N, 3] array
Unit normal vectors
'''
N = enforce_shape(N0, (3, ))
normal = []
for T in self.tangent:
N = vector.norm(N - T*vector.dot1(N, T))
normal.append(N)
normal = np.array(normal)
if self.closed and twist_to_close:
binormal = vector.cross(self.tangent, normal)
N = vector.norm(N - vector.dot1(N, self.tangent[0]))
angle = np.arctan2(vector.dot(binormal[0], N), vector.dot(normal[0], N))
phi = np.arange(len(self.path))*(-angle)/len(self.path)
normal = normal * np.cos(phi).reshape(-1, 1) + binormal * np.sin(phi).reshape(-1, 1)
return normal
[docs] def mesh_tube(self, width=None, outline=10):
'''
Create a Mesh object by tracing the Path, create a tube.
If the Path has a ``color`` attribute, the mesh will be colored
accordingly.
Uses :func:`ilpm.mesh.extrude_shape`.
Parameters
----------
width : [N] array, float or None, optional (default: 1)
The thickness of the tube/ribbon. If not specified, will be
given by the ``thickness`` attribute, or defaults to 1.
outline : int, optional (default: 10)
The number of points for the circumference of the
circle/ribbon tracing the outline.
'''
if width is None:
width = getattr(self, 'thickness', 1)
N0=np.eye(3)[np.argmin(self.path[0])]
N = self.parallel_transport_framing(N0=N0, twist_to_close=True)
B = vector.cross(self.tangent, N)
return mesh.extrude_shape(self.path, N, B, outline, closed=self.closed, scale=width)
[docs]def as_path(x):
'''Convert object to a path, if possible.'''
if isinstance(x, Path): return x
elif isinstance(x, (list, tuple, np.ndarray)): return Path(x)
elif isinstance(x, dict): return Path(**x)
else: raise ValueError("%s can't be converted to a Path" % x)
[docs]class Tangle(object):
'''Object which describes multiple 3D paths.
Most methods are the same as for :class:`Path`.
Any attribute which is valid for *all* of the child paths is also a valid
attribute of the Tangle.
(It will return a list of values from the child Paths.)
Parameters
----------
paths : list of Paths, or objects convertable to Paths
ignore_empty_paths : bool, optional (default: True)
If True, empty paths in a Tangle are ignored.
info : dict
Arbitrary info, saved in file
'''
def __init__(self, paths, ignore_empty_paths=True, info={}):
if isinstance(paths, (list, tuple)):
self.paths = map(as_path, paths)
elif isinstance(paths, Path):
self.paths = [paths]
if ignore_empty_paths:
self.paths = filter(lambda p: len(p), self.paths)
self.info = dict(info)
[docs] def copy(self):
'''Return a copt of the Tangle'''
return Tangle([p.copy() for p in self], info=self.info.copy())
[docs] def to_dict(self, path_omit=('derivative_points', 'refine_s')):
'''Returns a copy of the Tangle as a dict.
parameters
----------
path_omit : tuple of strings
Keywords to omit from the path sub-dictionaries.'''
return dict(
paths=map(lambda x: x.to_dict(omit=path_omit), self.paths),
info=self.info,
)
[docs] def save(self, fn):
'''Save the path as a file. The filename is not enforced, but ``.tangle``
is recommended.
The resulting file is a JSON version of ``self.to_dict``.
'''
d = self.to_dict()
d['__file_type__'] = 'JSON 3D TANGLE'
save_json(d, fn)
[docs] def compacted_segments(self):
'''Create a compacted path array and a number of points array.
Used primarily for OpenCL calculations.
Returns
-------
points : [total_segments + num_paths, 3] array
A ``vstack`` of all the paths in the tangle. If the path is closed,
the beginning point is added the end of the appropriate chunk and
``num_points = len(path) + 1``
num_segs : int array, shape [num_paths]
The number of segments (points) in each path.
path_i0, path_i1 : int array, shape [num_paths]
The beginning and end index of each path. There will be a gap of 1
between ``i1`` for path ``j`` and ``i0`` of path ``j+1`` because
``i1`` references the start point of the last segment.
(Exception to the rule: if a path is open and ``include_ends==True``,
then the last point will be included.)
'''
paths = [np.concatenate([p.path, p.path[:1]]) if p.closed else p.path for p in self.paths]
n_seg = np.array(map(len, paths))
ii = np.cumsum(n_seg)
return np.vstack(paths), n_seg-1, np.concatenate([[0], ii[:-1]]), ii-1
[docs] def compacted_segment_attribute(self, attr):
'''Return a compcated form of a attached attribute of the paths.
Designed to work in coordinate with :meth:`compacted_segments`.
Parameters
----------
attr : string
A valid attribute on all the paths (e.g. ``"tangent"``).
Returns
-------
a : array
A concatenated array of the attributes whose indices match those
returned by the paths from :meth:`compacted_segments`.
'''
a = [getattr(p, attr) for p in self.paths]
return np.concatenate([np.concatenate([a, a[:1]]) if p.closed else a for (a, p) in zip(a, self.paths)])
[docs] def rot(self, rv):
'''Rotate the tangle. (see :meth:`Path.rot`)'''
for path in self.paths: path.rot(rv)
[docs] def apply_basis(self, B):
'''Apply basis to the tangle. (see :meth:`Path.apply_basis`)'''
for path in self.paths: path.apply_basis(B)
[docs] def translate(self, offset):
'''Translate the tangle. (see :meth:`Path.translate`)'''
for path in self.paths: path.translate(offset)
[docs] def scale(self, scale):
'''Scale the tangle. (see :meth:`Path.scale`)'''
for path in self.paths: path.scale(scale)
def __getitem__(self, key):
'''Return one of the paths in the tangle.'''
return self.paths[key]
def __iter__(self):
'''Iterate over the paths in the tangle.'''
for path in self.paths:
yield path
def __len__(self):
'''Return the number of paths in the tangle.'''
return len(self.paths)
[docs] def decimate(self, ratio):
'''Decimate all the paths in the tangle. (see :meth:`Path.decimate`)'''
for path in self: path.decimate(ratio)
[docs] def total_points(self):
'''Return the total points in all paths'''
return sum(map(len, self))
[docs] def moment(self):
'''Return the first moment of the path. (see :meth:`Path.moment`)'''
return sum(map(Path.moment, self))
[docs] def orient(self, vec):
'''Orient tangle so specified vector is in the x-direction. (see :meth:`Path.orient`)'''
for path in self: path.orient(vec)
[docs] def center(self):
'''Return center of tangle. (see :meth:`Path.center`)'''
C = np.zeros(3)
w = 0
for path in self:
C += (path.seg_length * path.midpoint()).sum(0)
w += path.L
return C * (1./w)
[docs] def rms_radius(self):
'''Return the root mean squared radius of the tangle. (see :meth:`Path.rms_radius`)'''
R = 0
w = 0
C = self.center()
for path in self:
R += (path.seg_length * vector.mag1(path.midpoint() - C)**2).sum()
w += path.L
return np.sqrt(R / w)
[docs] def eliminate_close_points(self, min_dist):
'''Eliminate points which are too close. (see :meth:`Path.eliminate_close_points`)'''
if min_dist >= 0:
for path in self:
path.eliminate_close_points(min_dist)
[docs] def reduce_by_curvature(self, points_per_circle=100):
r'''Reduce a tangle by eliminating points in nearly-straight sections. (see :meth:`Path.reduce_by_curvature')'''
for p in self:
p.reduce_by_curvature(points_per_circle)
[docs] def twist(self):
'''Calculate the normal twist for each element in a tangle. (see :meth:`Path.twist`)
If you need to specify a normal vector other than the predefined
`normal`, call :meth:`Path.twist` on the elements of the Tangle.
'''
return map(Path.twist, self)
[docs] def twist_rate(self):
'''Calculate the normal twist rate for each element in a tangle. (see :meth:`Path.twist`)
If you need to specify a normal vector other than the predefined
`normal`, call :meth:`Path.twist` on the elements of the Tangle.
'''
return map(lambda t: Path.twist(t, per_point=True), self)
[docs] def total_length(self):
'''Compute the total length of all the paths in the Tangle.'''
return sum([path.L for path in self])
[docs] def crossing_matrix(self):
'''Compute the crossing number matrix for all paths in the Tangle.
Returns
-------
signed_crossings : [N, N] array
The total signed crossings for each pair of paths.
Self-crossing term is the Writhe, the off-diagonal terms are linking
numbers (integer-ness is enforced).
unsigned crossings : [N, N] array
The average number of unsigned crossings for each pair of paths,
averaged over all orientations.
'''
cpath, n_segs, i0, i1 = self.compacted_segments()
_check_code()
i0 = map(CL.int, i0)
i1 = map(CL.int, i1)
cl_path = CL.to_device(cpath.astype(CL.real))
local_mem = CL.local_memory(2 * int(np.dtype(CL.real).itemsize * CL.max_work_group_size))
scm = np.zeros((len(self), len(self)), CL.real)
ucm = np.zeros((len(self), len(self)), CL.real)
c = []
for j in range(len(self)):
#For path j we test crossings with all paths up to j
c.append(CL.zeros((i1[j], 2), CL.real))
#Queue up jobs
for i in range(j+1):
#__kernel void path_count_crossings(global const REAL* X,
# int i1, int i2, int j1, int j2,
# global REAL2 crossings, local REAL2 local_c)
CL.path_count_crossings(cl_path,
i0[i], i1[i], i0[j], i1[j],
c[-1], local_mem)
#Read data from jobs
for j in range(len(self)):
cj = c[j].get()
for i in range(j+1):
cc = cj[i0[i]:i1[i]].sum(0)
if i != j:
cc = np.round(cc)
scm[j, i] = cc[0]
ucm[j, i] = cc[1]
scm[i, j] = cc[0]
ucm[i, j] = cc[1]
return scm, ucm
[docs] def inductance_matrix(self, a=1.0, core_model=None, period=None):
r'''Compute the path inductance of a Tangle.
The path inductance is defined as:
:math:`\mathcal{E}_{ij} \approx \frac{1}{4 \pi} \oint \oint_{\left|\vec r_i - \vec r_j\right| > \frac{a}{2}} \frac{d\vec r_i \cdot d\vec r_j}{\left|\vec r - \vec r_j\right|}+ \delta_{ij} \frac{2-\alpha}{2\pi} L_i`
Parameters
----------
a : float, optional [default: 1.0]
The core size.
core_model : float or string, optional
The core model used to determine :math:`\alpha`; see :func:`set_default_core_model` for options.
If a float value is specified, this will be used to directly set :meth:`\alpha`.
If not specified, it is determined by the default set by :func:`set_default_core_model` [defaulut: ``"solid"``].
period : [3] array or float, optional
If specified, a 3x3x3 periodic copy of the original tangle is used to calculate inductances.
Returns
-------
inductance : (N, N) matrix
'''
if core_model is None:
alpha = CORE_MODELS[DEFAULT_CORE_MODEL]['alpha']
elif isinstance(core_model, basestring):
alpha = CORE_MODELS[core_model]['alpha']
else:
alpha = float(core_model)
Y = (2 - alpha) / (2*np.pi)
cpath, n_segs, i0, i1 = self.compacted_segments()
_check_code()
num_segs = len(cpath)
num_paths = len(i0)
if period is not None:
cp = []
ii0 = []
ii1 = []
num_copies = 0
for ox in [0, 1, -1]:
for oy in [0, 1, -1]:
for oz in [0, 1, -1]:
offset = np.array([ox, oy, oz], dtype='d') * period
# print offset
cp.append(cpath + offset)
# print cp[-1].shape, cp[-1].mean(0)
ii0.append(i0 + num_segs*num_copies)
ii1.append(i1 + num_segs*num_copies)
num_copies += 1
cpath = np.vstack(cp)
i0 = np.concatenate(ii0)
i1 = np.concatenate(ii1)
else:
num_copies = 1
i0 = map(CL.int, i0)
i1 = map(CL.int, i1)
cl_path = CL.to_device(cpath.astype(CL.real))
local_mem = CL.local_memory(int(np.dtype(CL.real).itemsize * CL.max_work_group_size))
cutoff = CL.real(0.5*a)
ind_pp = []
for j in range(len(self)):
#For path j we test crossings with all paths up to j
ind_pp.append(CL.zeros(i1[j], CL.real))
#Queue up jobs
for i in range(j+1):
# __kernel void path_inductance(global REAL* X,
# int i1, int i2, int j1, int j2, REAL cutoff,
# global REAL* inductance, local REAL* local_i)
for n in range(num_copies):
offset = n * num_paths
# print i0[i], i1[i], i0[j+offset], i1[j+offset]
CL.path_inductance(cl_path,
i0[i], i1[i], i0[j+offset], i1[j+offset], (cutoff if (i==j and n==0) else CL.real(0)),
ind_pp[-1], local_mem)
#Read data from jobs
ind = np.zeros((len(self), len(self)), CL.real)
for j in range(len(self)):
ind1 = ind_pp[j].get()
for i in range(j+1):
ind[i, j] = ind1[i0[i]:i1[i]].sum()
if i != j: ind[j, i] = ind[i, j]
else: ind[i, i] += self[i].L_polygon * Y
return ind
[docs] def closest_opposed(self, cutoff=0, per_point=False):
'''Find the closet point in the Tangle where paths are facing in opposite directions.
Tests the tangent vectors of all pairs of points in the Tangle, finding
the closest pair for which :math:``T_i \cdot T_j \leq cutoff``.
If no pairs meet this criterium, NaN is returned instead.
Parameters
----------
cutoff : float, optional (default: 0)
per_point : bool, optional (default: False)
If True, returns a list of arrays corresponding to the closest
distance from each individual point.
Returns
-------
r : float or array
The cloesest radius meeting the condition. (Optionally for each
individual point)
'''
_check_code()
Np = map(len, self.path)
Nt = sum(Np)
X = CL.to_device(np.concatenate(self.path).astype(CL.real))
T = CL.to_device(np.concatenate(self.tangent).astype(CL.real))
min_r = CL.to_device(-np.ones(Nt, dtype=CL.real))
local_mem = CL.local_memory(int(np.dtype(CL.real).itemsize * CL.max_work_group_size))
#__kernel void path_closest_opposed(
# global REAL* X, global REAL* T,
# int i1, int i2, int j1, int j2,
# REAL cutoff, char half_pair,
# global REAL* distance, local REAL* local_d)
CL.path_closest_opposed(
X, T,
0, Nt, 0, Nt,
CL.real(cutoff), CL.char(0 if per_point else 1),
min_r, local_mem
)
min_r = min_r.get()
min_r[np.where(min_r < 0)] = np.nan
if per_point:
ii = [0]+np.cumsum(Np).tolist()
return [min_r[ii[j]:ii[j+1]] for j in range(len(self))]
else:
return np.nanmin(min_r)
def __getattr__(self, k):
#This only gets called if an attribute is requested that has
# not already been computed.
vals = []
missing = 0
for p in self.paths:
if hasattr(p, k):
vals.append(getattr(p, k))
else:
missing += 1
if missing == len(self):
raise AttributeError("Tangle does not have property '%s'" % k)
elif missing:
raise AttributeError("Tangle does not have property '%s'\n(%d/%d of the child paths have it, but all are required)" % (k, missing, len(self)))
else:
return vals
[docs] def colorize(self, c, cmap='jet', clim=None):
'''Color a Tangle. (see :meth:`Path.colorize`)
The color may be specified either as a single item, or a list/tuple of
values. (e.g. a list of arrays that correspond to color values)
'''
if isinstance(c, basestring):
for p in self: p.colorize(c)
elif type(c) in (list, tuple):
if clim is None and not isintance(c[0], basestring):
c_all = np.concatenate(c)
clim = (np.nanmin(c_all), np.nanmax(c_all))
for cc, p in zip(c, self): p.colorize(cc, cmap=cmap, clim=clim)
else:
raise ValueError('color specification should either be single value or list of values')
[docs] def drop(self, attr):
'''Drop the specified attribute from the child paths.'''
for p in self:
if hasattr(p, attr): delattr(p, attr)
[docs] def mean(self, a):
'''Compute the mean of some variable over the tangle. See :meth:`Path.mean`'''
if isinstance(a, basestring):
a = [getattr(p, a) for p in self]
s = 0
L = 0
for aa, p in zip(a, self):
aa = enforce_shape(aa, (len(p), -1))
dL = p.seg_length_per_point()
s = s + (aa * dL).sum(0)
L += dL.sum()
return s / L
[docs] def color_by_proximity(self, p, c=None, w=None, rgb_only=True, power=-6, period=None):
'''Color a path by proximity to another path, tangle, or set of colored points. See :meth:`Path.color_by_proximity`'''
if hasattr(p, '_color_by_proximity_points'):
p, c, w = p._color_by_proximity_points()
p = enforce_shape(p, (-1, 3))
if c is None: raise ValueError('if specified as points, a color array must also be specified')
c = enforce_shape(c, (len(p), -1))
if c.shape[-1] == 3: c = np.concatenate([c, np.zeros((len(c), 1))], axis=-1)
elif c.shape[-1] != 4: raise ValueError('last dimension of color array must be 3 or 4')
if w is None: w = np.ones((len(c), 1))
else: w = enforce_shape(w, (len(c), 1) )
_check_code()
pw = CL.to_device(np.concatenate([p, w], axis=-1).astype(CL.real))
c = CL.to_device(c.astype(CL.real))
for path in self:
path._color_by_proximity_inner(pw, c, rgb_only, power, period=period)
def _color_by_proximity_points(self):
#Internal use
return np.concatenate(self.path), np.concatenate(self.color), np.concatenate([p.seg_length_per_point() for p in self])
[docs] def resample_with_spacing(self, dl, attached_interp='linear'):
'''Resample path evenly, spacing points by specified distance. See :meth:`Path.resample_with_spacing'
'''
paths = [p.resample_with_spacing(dl, attached_interp) for p in self]
return Tangle(paths, info=self.info)
[docs] def drop_shorter_than(self, N=None, L_polygon=None, L=None):
'''Remove short paths from the tangle, specified either in number of points and/or length.
Parameters
----------
N : int, optional
Clip by this number of points
L_polygon :
Clip by this length (from the ``L_polygon`` parameter of the path).
This avoids complete computation of the interpolations parameters,
which is much faster.
L : float, optional
Clip by this length (from the ``L`` parameter of the path)
'''
if N is not None:
self.paths = filter(lambda p: len(p.path) >= N, self.paths)
if L_polygon is not None:
self.paths = filter(lambda p: p.L_polygon >= L_polygon, self.paths)
if L is not None:
self.paths = filter(lambda p: p.L >= L, self.paths)
# def mesh_trace(self, width=None, outline=10, normal=None, ribbon_colors=None, ribbon_aspect=0.1):
# '''**NOT IMPLEMENTED YET**
#
# Create a Mesh object by tracing the Path, create a tube or ribbon.
#
# By default, creates a ribbon if the ``normal`` attribute is defined,
# and a tube if it is not.
#
# If the ``color`` attribute is present, the mesh will be colored
# accordingly.
#
# Uses :func:`ilpm.mesh.extrude_shape`.
#
# Parameters
# ----------
# width : [N] array, float or None, optional
# The thickness of the tube/ribbon. If not specified, will be
# given by the ``thickness`` attribute, or defaults to 1.
# outline : int or [M, 2] array, optional
# If int, the number of points for the circumference of the
# circle/ribbon tracing the outline.
# If array, the outline shape (first axis points in the normal
# direction).
# normal : [N, 3] array or False, optional
# The normal vector to use for the tracing. Note that this normal
# vector will be passed through :meth:`gen_normal`; if non-orthogonal
# normal vectors are required, use :func:`ilpm.mesh.extrude_shape`
# directly.
# If the normal is not specified, it will be derived from the
# ``normal`` field, or if not present a circular tube will be created
# instead.
# To force a circular tube when the ``normal`` is present, use
# ``normal==False``
# ribbon_colors : list of 2 color specifications, optional
# If specified, and a ribbon is being created, these colors will
# appear on the front and back of the ribbon. Ignored for tubes.
# ribbon_aspect : float, optional (default: 0.1)
# The thickness of the ribbon with respect to the width.
# Ignored for tubes.
# '''
#
# raise ImportError('Not implemented yet -- sorry!')
#
# if width is None:
# width = getattr(self, 'thickness', 1)
#
# if normal is None:
# if hasattr('normal', self):
# pass
[docs] def mesh_tube(self, width=None, outline=10):
'''
Create a Mesh object by tracing the Tangle, create a tube.
See :meth:`Path.mesh_tube` for details.
'''
m = None
for path in self:
mm = path.mesh_tube(width=width, outline=outline)
if m is None: m = mm
else: m += mm
return m
[docs] def reduce_topology(self):
'''Remove as many points from the tangle as possible without changing the crossings.
The remaining tangle will have the same topology, but with a small
number of points.
Topologically trivial loops will be removed; if the entire tangle is
topologically trivial it should create an empty tangle.
The algorithm works by choosing a random point and testing if it can
be removed without changing any crossings. (It tests for this by
looking for intersections between segments and the triangle formed by
the three-point neighborhood of the point.)
*Note:* It is possible that the reduction may get stuck if the tangled
is very tangled. Typically, however, it is able to reduce the tangle
succesfully.
*Note 2:* This algorithm can be slow; it is highly recommended to run
:meth:`reduce_by_curvature` on the tangle first.
'''
for tries in range(self.total_points()):
if not len(self):
#print "It's trivial!"
break
#Choose a random point to start the reduction
npp = map(len, self)
n_tot = sum(npp)
for i in ((np.arange(n_tot) + np.random.randint(0, n_tot)) % n_tot):
#Find which loop this is in an which point it corresponds to
ip = 0
while i >= npp[ip]:
i -= npp[ip]
ip += 1
#Three points form a triangle: we will now test if any other points
# pierce this triangle; if not, the mid point can be removed without
# changing the topology.
Xt = self[ip].path[(i+np.arange(-1, 2)) % len(self[ip])]
#Form an orthnormal basis with the triangle normal as X and the first
# vector as Y, then subtract the basis point
B = vector.normalize_basis([vector.cross(Xt[2] - Xt[1], Xt[1] - Xt[0]), Xt[2] - Xt[1]])
Xt = vector.apply_basis(Xt, B)
X0 = Xt[1].copy()
Xt -= X0
#Find all segments which cross the triangle plane
P = []
for jp, path in enumerate(self.path):
if ip == jp:
#If this is the same path, there are _four_ segments that can't
# posibly cross the triangle!
S0 = np.roll(path, 2-i, axis=0)[4:]
S1 = np.roll(path, 1-i, axis=0)[4:]
else:
S0 = path
S1 = np.roll(path, -1, axis=0)
#Get X coordinate of two points on segment
x0 = vector.dot(S0, B[0]) - X0[0]
x1 = vector.dot(S1, B[0]) - X0[0]
#First test gets the edge cases, but prevents zero divisions;
# completely in plane segments are just ignored, which is the only
# sensible course.
crosses_plane = np.where(((x0 == 0) & (x1 != 0)) | (x0 * x1 < 0))[0]
if len(crosses_plane):
x0 = x0[crosses_plane]
x1 = x1[crosses_plane]
#Compute point in plane crossing
s = ((0 - x0) / (x1 - x0))[:, np.newaxis]
Xp = S0[crosses_plane]*(1-s) + S1[crosses_plane]*s
P.append((vector.apply_basis(Xp, B) - X0)[:, 1:])
if P:
x, y = np.concatenate(P).T
#Convert to barycentric coordinates
# http://en.wikipedia.org/wiki/Barycentric_coordinate_system
#2D points of triangle
((x1, y1), (x2, y2), (x3, y3)) = Xt[:, 1:]
div = 1./( (y2-y3)*(x1-x3) + (x3-x2)*(y1-y3) )
if not np.isfinite(div):
warnings.warn('encountered null triangle in Tangle.reduce_topology; this may be fixable by running Tangle.reduce_by_curvature first.', ShapeWarning, stacklevel=2)
continue
a = ((y2-y3)*(x-x3) + (x3-x2)*(y-y3)) * div
b = ((y3-y1)*(x-x3) + (x1-x3)*(y-y3)) * div
no_crosses = ((a < 0) | (b < 0) | ((a+b) > 1)).all()
else:
no_crosses = True
#We get to remove a point!
if no_crosses:
if len(self[ip]) <= 3: #This path is done!
self.paths.pop(ip)
else:
self[ip].keep(np.concatenate([np.arange(0, i), np.arange(i+1, len(self[ip]))]))
break #Restart the random search
else: #We searched every point, to no avail!
#print "It's NOT trivial!"
#self.save(os.path.join(od, 'pr_%08d.tangle'%(tries+1)))
break
class NumpyAwareJSONEncoder(json.JSONEncoder):
def default(self, obj):
if isinstance(obj, np.ndarray):
return obj.tolist()
return json.JSONEncoder.default(self, obj)
def save_json(obj, fn):
with open(fn, 'wt') as f:
json.dump(obj, f, cls=NumpyAwareJSONEncoder, sort_keys=True, indent=2, separators=(',', ': '))
class FileTypeError(Exception):
pass
[docs]def load_path(fn, ignore_empty_paths=True):
'''Load a path or tangle from a file.
Supports either the new style path/tangle format as well as the older
".json" tangle format.
Parameters
----------
fn : string
A filename corresponding to a valid path.
ignore_empty_paths : bool, optional (default: True)
If True, empty paths in a Tangle are ignored.
Returns
-------
path : Path or Tangle
'''
with open(fn, 'rt') as f:
dat = json.load(f)
ft = dat.pop('__file_type__', None)
if ft == 'JSON 3D PATH':
return Path(**dat)
elif ft == 'JSON 3D TANGLE':
return Tangle(ignore_empty_paths=ignore_empty_paths, **dat)
else:
if 'traces' in dat:
paths = dat.pop('traces')
return Tangle(paths, info=dat, ignore_empty_paths=ignore_empty_paths)
else:
raise FileTypeError('%s does not appear to be a 3D path or tangle file' % fn)
[docs]def load_tangle(fn, ignore_empty_paths=True):
'''Similiar :meth:`load_path`, but always returns a tangle object.
Parameters
----------
fn : string
A filename corresponding to a valid path or tangle.
ignore_empty_paths : bool, optional (default: True)
If True, empty paths in a Tangle are ignored.
Returns
-------
tangle : Tangle
If the filename corresponds to a Path, it is converted to a Tangle with
a single path.
'''
T = load_path(fn, ignore_empty_paths=ignore_empty_paths)
if isinstance(T, Path):
T = Tangle([Path])
return T
[docs]def torus_knot(p, q, r=1, a=1./3, NP=128):
'''Generate a (generalized) torus knot.
Constructed so that ``p`` is the number of filaments in a bundle
cross-section and ``p/q`` is the twist number per revolution.
Parameters
----------
p, q: integer
Winding numbers
r, a: float, optional
Major and minor radius (default: r=1, a=1/3)
NP: int, optional
The number of points per revolution (default: 128)
Returns
-------
tangle : Tangle
A tangle of the torus knot.
'''
Nf = fractions.gcd(p, q)
p0 = p / Nf;
phi = np.arange(NP * p0) * 2*np.pi / NP
c = np.cos(phi)
s = np.sin(phi)
paths = []
for n in range(Nf):
theta = phi * q/p + 2*np.pi * n/p
rr = r + a * np.cos(theta)
paths.append(np.array([rr * c, rr * s, a * np.sin(theta)]).T)
return Tangle(paths, info={'description':'Torus knot, r_0=%f, a=%f, p=%d, q=%d' % (r, a, p, q)})
if __name__ == '__main__':
N = 100
closed = True
if closed:
phi = np.arange(N) * 2*np.pi / N
else:
phi = np.arange(N+1) * np.pi / N
#phi += 0.2 * np.cos(phi)
x = vector.vec(x=np.cos(phi), y=np.sin(phi))
test = Path(x, closed=closed, derivative_points=5)
print test.moment()
print
#print isinstance(test, Path)
#print type(x), isinstance(x, (list, tuple, np.ndarray))
test_t = Tangle([test.path.copy(), x+5])
print test_t.moment()
print test_t.center()
print test_t.rms_radius()
#test_t.rot((0, np.pi/2, 0))
test_t.orient(test_t.moment())
print 'oriented by moment'
print test_t.moment()
print test_t.center()
print test_t.rms_radius()
rv = np.array([0, np.pi/2, 0])
print vector.apply_basis(np.eye(3), vector.rot(rv))
print vector.rot(rv, np.eye(3))
X, Y, Z = np.eye(3)
print vector.apply_basis(np.eye(3), vector.normalize_basis([Z]))
sys.exit()
print (phi[1] - phi[0]) / test.s[1]
# print test.s.shape
# print test.L
# print test._Ds(-1)
# print test.ds_weights
print vector.mag(test.partial_s('path'))
#print test._Ds(1)
test.save('test.path')
test_t.save('test.tangle')
test2 = load_path('test.path')
test2_t = load_path('test.tangle')
print (test.path == test.path).all()
print type(test2_t)
