From e0211e21f3db0bb27a04490d712dabe04c0fed26 Mon Sep 17 00:00:00 2001 From: mdroe Date: Fri, 13 Jun 2014 22:04:04 +0000 Subject: More speedups. git-svn-id: http://svn.stsci.edu/svn/ssb/stsci_python/stsci.sphere/trunk@32541 fe389314-cf27-0410-b35b-8c050e845b92 Former-commit-id: 36e70c1cac043d59089f386a6b6d730c4887fa84 --- src/math_util.c | 90 +++++++++++++++++++++++++++++------------ stsci/sphere/graph.py | 14 +++++-- stsci/sphere/polygon.py | 31 +++++++------- stsci/sphere/test/test.py | 5 ++- stsci/sphere/test/test_union.py | 2 +- stsci/sphere/vector.py | 62 ++++++++++++++-------------- 6 files changed, 124 insertions(+), 80 deletions(-) diff --git a/src/math_util.c b/src/math_util.c index 9dbb32a..2aae11a 100644 --- a/src/math_util.c +++ b/src/math_util.c @@ -49,6 +49,9 @@ typedef npy_intp intp; INIT_OUTER_LOOP_4 \ intp s4 = *steps++; +#define BEGIN_OUTER_LOOP_2 \ + for (N_ = 0; N_ < dN; N_++, args[0] += s0, args[1] += s1) { + #define BEGIN_OUTER_LOOP_3 \ for (N_ = 0; N_ < dN; N_++, args[0] += s0, args[1] += s1, args[2] += s2) { @@ -96,13 +99,6 @@ save_pointl(const long double* in, char *out, const intp s) { *(double *)out = (double)in[2]; } -static inline void -cross(const double *A, const double *B, double *C) { - C[0] = A[1]*B[2] - A[2]*B[1]; - C[1] = A[2]*B[0] - A[0]*B[2]; - C[2] = A[0]*B[1] - A[1]*B[0]; -} - static inline void crossl(const long double *A, const long double *B, long double *C) { C[0] = A[1]*B[2] - A[2]*B[1]; @@ -111,13 +107,17 @@ crossl(const long double *A, const long double *B, long double *C) { } static inline void -normalize(double *A) { +normalize_output(long double *A, double *B) { double l = A[0]*A[0] + A[1]*A[1] + A[2]*A[2]; if (l != 1.0) { l = sqrt(l); - A[0] /= l; - A[1] /= l; - A[2] /= l; + B[0] = A[0] / l; + B[1] = A[1] / l; + B[2] = A[2] / l; + } else { + B[0] = A[0]; + B[1] = A[1]; + B[2] = A[2]; } } @@ -206,6 +206,35 @@ static PyUFuncGenericFunction inner1d_functions[] = { DOUBLE_inner1d }; static void * inner1d_data[] = { (void *)NULL }; static char inner1d_signatures[] = { PyArray_DOUBLE, PyArray_DOUBLE, PyArray_DOUBLE }; +/*/////////////////////////////////////////////////////////////////////////// + normalize +*/ + +char *normalize_signature = "(i)->(i)"; + +static void +DOUBLE_normalize(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)) +{ + long double IN[3]; + double OUT[3]; + + INIT_OUTER_LOOP_2 + intp is1=steps[0], is2=steps[1]; + BEGIN_OUTER_LOOP_2 + char *ip1=args[0], *op=args[1]; + + load_pointl(ip1, is1, IN); + + normalize_output(IN, OUT); + + save_point(OUT, op, is2); + END_OUTER_LOOP +} + +static PyUFuncGenericFunction normalize_functions[] = { DOUBLE_normalize }; +static void * normalize_data[] = { (void *)NULL }; +static char normalize_signatures[] = { PyArray_DOUBLE, PyArray_DOUBLE }; + /*/////////////////////////////////////////////////////////////////////////// cross */ @@ -215,21 +244,21 @@ char *cross_signature = "(i),(i)->(i)"; static void DOUBLE_cross(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)) { - double A[3]; - double B[3]; - double C[3]; + long double A[3]; + long double B[3]; + long double C[3]; INIT_OUTER_LOOP_3 intp is1=steps[0], is2=steps[1], is3=steps[2]; BEGIN_OUTER_LOOP_3 char *ip1=args[0], *ip2=args[1], *op=args[2]; - load_point(ip1, is1, A); - load_point(ip2, is2, B); + load_pointl(ip1, is1, A); + load_pointl(ip2, is2, B); - cross(A, B, C); + crossl(A, B, C); - save_point(C, op, is3); + save_pointl(C, op, is3); END_OUTER_LOOP } @@ -250,22 +279,22 @@ char *cross_and_norm_signature = "(i),(i)->(i)"; static void DOUBLE_cross_and_norm(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)) { - double A[3]; - double B[3]; - double C[3]; + long double A[3]; + long double B[3]; + long double C[3]; INIT_OUTER_LOOP_3 intp is1=steps[0], is2=steps[1], is3=steps[2]; BEGIN_OUTER_LOOP_3 char *ip1=args[0], *ip2=args[1], *op=args[2]; - load_point(ip1, is1, A); - load_point(ip2, is2, B); + load_pointl(ip1, is1, A); + load_pointl(ip2, is2, B); - cross(A, B, C); - normalize(C); + crossl(A, B, C); + normalizel(C); - save_point(C, op, is3); + save_pointl(C, op, is3); END_OUTER_LOOP } @@ -454,6 +483,15 @@ addUfuncs(PyObject *dictionary) { PyDict_SetItemString(dictionary, "inner1d", f); Py_DECREF(f); + f = PyUFunc_FromFuncAndDataAndSignature( + normalize_functions, normalize_data, normalize_signatures, 2, 1, 1, + PyUFunc_None, "normalize", + "inner on the last dimension and broadcast on the rest \n" \ + " \"(i)->(i)\" \n", + 0, normalize_signature); + PyDict_SetItemString(dictionary, "normalize", f); + Py_DECREF(f); + f = PyUFunc_FromFuncAndDataAndSignature( cross_functions, cross_data, cross_signatures, 2, 2, 1, PyUFunc_None, "cross", diff --git a/stsci/sphere/graph.py b/stsci/sphere/graph.py index 0fbc9a7..bae1c60 100644 --- a/stsci/sphere/graph.py +++ b/stsci/sphere/graph.py @@ -270,11 +270,17 @@ class Graph: # Don't add nodes that already exist. Update the existing # node's source_polygons list to include the new polygon. - point = vector.normalize_vector(*point) + point = vector.normalize_vector(point) - # TODO: Vectorize this - for node in self._nodes: - if np.all(np.abs(node._point - point) < 2 ** -32): + if len(self._nodes): + nodes = list(self._nodes) + node_array = np.array([node._point for node in nodes]) + + diff = np.all(np.abs(point - node_array) < 2 ** -32, axis=-1) + + indices = np.nonzero(diff)[0] + if len(indices): + node = nodes[indices[0]] node._source_polygons.update(source_polygons) return node diff --git a/stsci/sphere/polygon.py b/stsci/sphere/polygon.py index f7c44be..5167171 100644 --- a/stsci/sphere/polygon.py +++ b/stsci/sphere/polygon.py @@ -168,11 +168,11 @@ class SphericalPolygon(object): # Convert to Cartesian x, y, z = vector.radec_to_vector(ra, dec, degrees=degrees) + points = np.dstack((x, y, z))[0] + if center is None: - xc = x.mean() - yc = y.mean() - zc = z.mean() - center = vector.normalize_vector(xc, yc, zc) + center = points.mean(axis=0) + vector.normalize_vector(center, output=center) else: center = vector.radec_to_vector(*center, degrees=degrees) @@ -556,24 +556,21 @@ class SphericalPolygon(object): # Rotate polygon so that center of polygon is at north pole centroid = np.mean(points[:-1], axis=0) - centroid = vector.normalize_vector(*centroid) + centroid = vector.normalize_vector(centroid) points = self._points - (centroid + np.array([0, 0, 1])) - vector.normalize_vector( - points[:, 0], points[:, 1], points[:, 2], inplace=True) + vector.normalize_vector(points, output=points) - X = [] - Y = [] + XYs = [] for A, B in zip(points[:-1], points[1:]): length = great_circle_arc.length(A, B, degrees=True) interp = great_circle_arc.interpolate(A, B, length * 4) - x, y, z = vector.normalize_vector( - interp[:, 0], interp[:, 1], interp[:, 2], inplace=True) - x, y = vector.equal_area_proj(x, y, z) - X.extend(x) - Y.extend(y) - - X = np.array(X) - Y = np.array(Y) + vector.normalize_vector(interp, output=interp) + XY = vector.equal_area_proj(interp) + XYs.append(XY) + + XY = np.vstack(XYs) + X = XY[..., 0] + Y = XY[..., 1] return np.abs(np.sum(X[:-1] * Y[1:] - X[1:] * Y[:-1]) * 0.5 * np.pi) diff --git a/stsci/sphere/test/test.py b/stsci/sphere/test/test.py index 8e6e270..5ab7bbb 100644 --- a/stsci/sphere/test/test.py +++ b/stsci/sphere/test/test.py @@ -19,8 +19,9 @@ graph.DEBUG = True def test_normalize_vector(): x, y, z = np.ogrid[-100:100:11,-100:100:11,-100:100:11] - xn, yn, zn = vector.normalize_vector(x.flatten(), y.flatten(), z.flatten()) - l = np.sqrt(xn ** 2 + yn ** 2 + zn ** 2) + xyz = np.dstack((x.flatten(), y.flatten(), z.flatten()))[0] + xyzn = vector.normalize_vector(xyz) + l = np.sqrt(np.sum(xyzn * xyzn, axis=-1)) assert_almost_equal(l, 1.0) diff --git a/stsci/sphere/test/test_union.py b/stsci/sphere/test/test_union.py index 60ea572..a6e5793 100644 --- a/stsci/sphere/test/test_union.py +++ b/stsci/sphere/test/test_union.py @@ -209,7 +209,7 @@ def test_difficult_unions(): poly = polygon.SphericalPolygon(points, inside) polys.append(poly) - polygon.SphericalPolygon.multi_union(polys) + polygon.SphericalPolygon.multi_union(polys[:len(polys)/2]) if __name__ == '__main__': diff --git a/stsci/sphere/vector.py b/stsci/sphere/vector.py index c6887c4..bb444e6 100644 --- a/stsci/sphere/vector.py +++ b/stsci/sphere/vector.py @@ -42,6 +42,12 @@ from __future__ import unicode_literals # THIRD-PARTY import numpy as np +try: + from . import math_util + HAS_C_UFUNCS = True +except ImportError: + HAS_C_UFUNCS = False + __all__ = ['radec_to_vector', 'vector_to_radec', 'normalize_vector', 'rotate_around'] @@ -121,9 +127,9 @@ def vector_to_radec(x, y, z, degrees=True): \delta = \arctan2(z, \sqrt{x^2 + y^2}) """ - x = np.asanyarray(x) - y = np.asanyarray(y) - z = np.asanyarray(z) + x = np.asanyarray(x, dtype=np.float64) + y = np.asanyarray(y, dtype=np.float64) + z = np.asanyarray(z, dtype=np.float64) result = ( np.arctan2(y, x), @@ -135,39 +141,37 @@ def vector_to_radec(x, y, z, degrees=True): return result -def normalize_vector(x, y, z, inplace=False): +def normalize_vector(xyz, output=None): r""" Normalizes a vector so it falls on the unit sphere. - *x*, *y*, *z* may be scalars or 1-D arrays - Parameters ---------- - x, y, z : scalars or 1-D arrays of the same length + xyz : Nx3 array of vectors The input vectors - inplace : bool, optional - When `True`, the original arrays will be normalized in place, - otherwise a normalized copy is returned. + output : Nx3 array of vectors, optional + The array to store the results in. If `None`, a new array + will be created and returned. Returns ------- - X, Y, Z : scalars of 1-D arrays of the same length - The normalized output vectors + output : Nx3 array of vectors """ - x = np.asanyarray(x) - y = np.asanyarray(y) - z = np.asanyarray(z) + xyz = np.asanyarray(xyz, dtype=np.float64) - l = np.sqrt(x ** 2 + y ** 2 + z ** 2) + if output is None: + output = np.empty(xyz.shape, dtype=np.float64) - if inplace: - x /= l - y /= l - z /= l - return (x, y, z) - else: - return (x / l, y / l, z / l) + if HAS_C_UFUNCS: + math_util.normalize(xyz, output) + return output + + l = np.sqrt(np.sum(xyz * xyz, axis=-1)) + + output = xyz / np.expand_dims(l, 2) + + return output def rotate_around(x, y, z, u, v, w, theta, degrees=True): @@ -212,19 +216,19 @@ def rotate_around(x, y, z, u, v, w, theta, degrees=True): return X, Y, Z -def equal_area_proj(x, y, z): +def equal_area_proj(points): """ Transform the coordinates to a 2-dimensional plane using the Lambert azimuthal equal-area projection. Parameters ---------- - x, y, z : scalars or 1-D arrays + points : Nx3 array of vectors The input vectors Returns ------- - X, Y : tuple of scalars or arrays of the same length + output : Nx2 array of points Notes ----- @@ -237,7 +241,5 @@ def equal_area_proj(x, y, z): Y = \sqrt{\frac{2}{1-z}}y """ - scale = np.sqrt(2.0 / (1.0 - z)) - X = scale * x - Y = scale * y - return X, Y + scale = np.sqrt(2.0 / (1.0 - points[..., 2])) + return np.expand_dims(scale, 2) * points[:, 0:2] -- cgit