Initial move over from core (#1)

There's a bunch of general math utilities in core that Mesh will depend on, core depends on. More, there's a bunch of stuff that depends on core in part just because it needs these generic math utilities. So this seemed like a good first thing to break out.

This PR is really all structural, there are no interesting changes in the actual implementations.

Author: algrs

RENAME-TO-LICENSE-IF-PUBLIC-ELSE-DELETE LICENSERENAME-TO-LICENSE-IF-PUBLIC-ELSE-DELETE renamed to LICENSE

@@ -1,75 +1,122 @@
-skeleton-py
-===========
-
-A project starter for Body Labs Python projects.
-
-To use this skeleton:
-
-1. Clone the skeleton:
-
-        git clone [email protected]:bodylabs/template-py.git my-library
-
-2. Re-initialize the git repository:
-
-        cd my-library && rm -rf .git && git init
-
-3. Publish the license, or delete it, depending whether or not the project is
-   open source:
-
-        mv RENAME-TO-LICENSE-IF-PUBLIC-ELSE-DELETE LICENSE
-        rm RENAME-TO-LICENSE-IF-PUBLIC-ELSE-DELETE
-
-4. Commit the changed files.
-
-        git add . && git commit -m "Project skeleton"
-
-5. Make updates as needed.
-
-
-- - - - - - - - - - - - -
-
-example
-=======
-
-Short description of the module.
-
-
-Features
---------
-
-- This
-- That
-- The other
-
-
-Examples
---------
-
-```py
-from example.hello import hello
-hello()
-```
-
-```sh
-hello
-```
+blmath
+======
+
+A collection of math related utilities used by many bits of BodyLabs code.
+
+Requirements
+------------
+On macOS:
+
+    brew install homebrew/science/suite-sparse
+    brew install homebrew/science/opencv --without-numpy
+
+On Linux:
+
+    sudo apt-get install python-opencv libsuitesparse-dev
+
+On windows:
+
+TODO: Windows install nstructions
+
+
+blmath.numerics
+---------------
+
+Functions for manipulating numeric arrays, numbers, and linear algebra.
+
+The [most commonly used of these](__init__.py) are directly imported into
+`blmath.numerics`.
+
+- [blmath.numerics.vx](vector_shortcuts.py) is a namespace of common linear
+  algebra operations. These are easily expressed in numpy, but abstracted for
+  readability purposes.
+- [blmath.numerics.coercion](coercion.py) contains a validation function
+  `as_numeric_array`, which produces useful error messages up front on bad
+  inputs, in place of cryptic messages like "cannot broadcast..." later on.
+- [blmath.numerics.operations](operations.py) contains basic numerical
+  operations such as `zero_safe_divide`.
+- [blmath.numerics.predicates](predicates.py) contains functions like
+  `isnumeric`.
+- [blmath.numerics.rounding](rounding.py) contains functions including
+  "round to nearest" and `roundedlist`.
+- [blmath.numerics.numpy_ext](numpy_ext.py) contains numpy utility
+  functions.
+- [blmath.numerics.matlab](matlab.py) contains some matlab shortcuts which
+  have no numpy equivalent. At MPI the fitting code was originally written in
+  Matlab before it was ported to Python.
+
+[blmath.numerics.linalg](linalg) contains linear algebra operations.
+
+- [blmath.numerics.linalg.sparse_cg](linalg/sparse_cg.py) contains a faster
+  matrix solve optimized for sparse Jacobians.
+- [blmath.numerics.linalg.lchol](linalg/lchol.py) contains a Cythonized
+  implementation of Cholesky factorization.
+- [blmath.numerics.linalg.isomorphism](linalg/isomorphism.py) computes the
+  isomorphism between two bases.
+- [blmath.numerics.linalg.gram_schmidt](linalg/gram_schmidt.py) provides a
+  function for orthonormalization.
+
+blmath.geometry
+---------------
+
+Geometric operations, transforms, and primitives, in 2D and 3D.
+
+The [most commonly used of these](__init__.py) are directly imported into
+`blmath.geometry`.
+
+- [blmath.geometry.Box](primitives/box.py) represents an axis-aligned
+  cuboid.
+- [blmath.geometry.Plane](primitives/plane.py) represents a 2-D plane in
+  3-space (not a hyperplane).
+- [blmath.geometry.Polyline](primitives/polyline.py) represents an
+  unconstrained polygonal chain in 3-space.
+
+`blmath.geometry.transform` includes code for 3D transforms.
+
+- [blmath.geometry.transform.CompositeTransform](transform/composite.py)
+  represents a composite transform using homogeneous coordinates. (Thanks avd!)
+- [blmath.geometry.transform.CoordinateManager](transform/coordinate_manager.py)
+  provides a convenient interface for named reference frames within a stack of
+  transforms and projecting points from one reference frame to another.
+- [blmath.geometry.transform.find_rigid_transform](transform/rigid_transform.py)
+  finds a rotation and translation that closely transforms one set of points to
+  another. Its cousin `find_rigid_rotation` does the same, but only allows
+  rotation, not translation.
+- [blmath.geometry.transform.rotation.rotation_from_up_and_look](transform/rotation.py)
+  produces a rotation matrix that gets a mesh into the canonical reference frame
+  from "up" and "look" vectors.
+
+Other modules:
+
+- [blmath.geometry.apex](apex.py) provides functions for finding the most
+  extreme point.
+- [blmath.geometry.barycentric](barycentric.py) provides a function for
+  projecting a point to a triangle using barycentric coordinates.
+- [blmath.geometry.convexify](convexify.py) provides a function for
+  producing a convex hull from a mostly-planar curve.
+- [blmath.geometry.segment](segment.py) provides functions for working with
+  line segments in n-space.
+
+blmath.units
+------------
+TODO write something here
 
 
 Development
 -----------
 
 ```sh
 pip install -r requirements_dev.txt
-rake test
+rake unittest
 rake lint
 ```
 
 
 Contribute
 ----------
 
-- Issue Tracker: github.com/bodylabs/example/issues
-- Source Code: github.com/bodylabs/example
+- Issue Tracker: github.com/bodylabs/blmath/issues
+- Source Code: github.com/bodylabs/blmath
 
 Pull requests welcome!
 

README.md

@@ -41,7 +41,7 @@ task :unittest do
 end
 
 task :lint => :require_style_config do
-  raise unless system "bodylabs-python-style/bin/pylint_test example --min_rating 10.0"
+  raise unless system "bodylabs-python-style/bin/pylint_test blmath --min_rating 10.0"
 end
 
 desc "Remove .pyc files"

Rakefile

@@ -0,0 +1,18 @@
+import os
+from baiji.pod import AssetCache, VersionedCache
+from baiji.pod.config import Config
+
+config = Config()
+config.CACHE_DIR = os.path.expanduser('~/.bodylabs_static_cache')
+config.IMMUTABLE_BUCKETS = ['bodylabs-versioned-assets', 'bodylabs-versioned-assets-tokyo']
+config.DEFAULT_BUCKET = os.getenv('ASSET_BUCKET', 'bodylabs-assets')
+
+sc = AssetCache(config)
+
+MANIFEST_PATH = os.path.join(os.path.dirname(__file__), '..', 'manifest.json')
+BUCKET = os.getenv('VC_BUCKET', 'bodylabs-versioned-assets')
+
+vc = VersionedCache(
+    cache=sc,
+    manifest_path=MANIFEST_PATH,
+    bucket=BUCKET)

bin/hello

@@ -0,0 +1,3 @@
+from blmath.geometry.primitives.box import Box
+from blmath.geometry.primitives.plane import Plane
+from blmath.geometry.primitives.polyline import Polyline

example/__init__.py blmath/__init__.pyexample/__init__.py renamed to blmath/__init__.py

@@ -0,0 +1,57 @@
+import numpy as np
+from blmath.numerics import vx
+
+def apex(points, axis):
+    '''
+    Find the most extreme point in the direction of the axis provided.
+
+    axis: A vector, which is an 3x1 np.array.
+
+    '''
+    coords_on_axis = points.dot(axis)
+    return points[np.argmax(coords_on_axis)]
+
+def inflection_points(points, axis, span):
+    '''
+    Find the list of vertices that preceed inflection points in a curve. The curve is differentiated
+    with respect to the coordinate system defined by axis and span.
+
+    axis: A vector representing the vertical axis of the coordinate system.
+    span: A vector representing the the horiztonal axis of the coordinate system.
+
+    returns: a list of points in space corresponding to the vertices that
+    immediately preceed inflection points in the curve
+    '''
+
+    coords_on_span = points.dot(span)
+    dx = np.gradient(coords_on_span)
+    coords_on_axis = points.dot(axis)
+
+    # Take the second order finite difference of the curve with respect to the
+    # defined coordinate system
+    finite_difference_2 = np.gradient(np.gradient(coords_on_axis, dx), dx)
+
+    # Compare the product of all neighboring pairs of points in the second derivative
+    # If a pair of points has a negative product, then the second derivative changes sign
+    # at one of those points, signalling an inflection point
+    is_inflection_point = [finite_difference_2[i] * finite_difference_2[i + 1] <= 0 for i in range(len(finite_difference_2) - 1)]
+
+    inflection_point_indices = [i for i, b in enumerate(is_inflection_point) if b]
+
+    if len(inflection_point_indices) == 0:
+        return []
+
+    return points[inflection_point_indices]
+
+def farthest(from_point, to_points):
+    '''
+    Find the farthest point among the inputs, to the given point.
+
+    Return a tuple: farthest_point, index_of_farthest_point.
+    '''
+    absolute_distances = vx.magnitude(to_points - from_point)
+
+    index_of_farthest_point = np.argmax(absolute_distances)
+    farthest_point = to_points[index_of_farthest_point]
+
+    return farthest_point, index_of_farthest_point

blmath/cache.py

@@ -0,0 +1,43 @@
+import numpy as np
+
+def barycentric_coordinates_of_projection(p, q, u, v):
+    """ Given a point, gives projected coords of that point to a triangle
+    in barycentric coordinates.
+
+    See Heidrich, Computing the Barycentric Coordinates of a Projected Point, JGT 05
+    at http://www.cs.ubc.ca/~heidrich/Papers/JGT.05.pdf
+
+    Args:
+        p: point to project
+        q: a vertex of the triangle to project into
+        u,v: edges of the the triangle such that it has vertices q, q+u, q+v
+
+    Returns:
+        b: barycentric coordinates of p's projection in triangle defined by q,u,v
+            vectorized so p,q,u,v can all be 3xN
+    """
+
+    p = p.T
+    q = q.T
+    u = u.T
+    v = v.T
+
+    n = np.cross(u, v, axis=0)
+    s = np.sum(n*n, axis=0)
+
+    # If the triangle edges are collinear, cross-product is zero,
+    # which makes "s" 0, which gives us divide by zero. So we
+    # make the arbitrary choice to set s to epsv (=numpy.spacing(1)),
+    # the closest thing to zero
+    if np.isscalar(s):
+        s = s if s else np.spacing(1)
+    else:
+        s[s == 0] = np.spacing(1)
+
+    oneOver4ASquared = 1.0 / s
+    w = p - q
+    b2 = np.sum(np.cross(u, w, axis=0) * n, axis=0) * oneOver4ASquared
+    b1 = np.sum(np.cross(w, v, axis=0) * n, axis=0) * oneOver4ASquared
+    b = np.vstack((1 - b1 - b2, b1, b2))
+
+    return b.T

blmath/geometry/__init__.py

@@ -0,0 +1,47 @@
+def convexify_planar_curve(polyline, flatten=False, want_vertices=False):
+    '''
+    Take the convex hull of an almost planar 2-D curve in three-space.
+
+    Params:
+        - polyline:
+            An instance of Polyline.
+        - flatten:
+            Boolean; if True, rotated curve will be flattened
+            on the y-axis, thereby losing length. This loss
+            may not be offset by the convex hull.
+        - want_vertices:
+            Boolean; do you want the indices to the convex hull vertices?
+    '''
+    import numpy as np
+    from scipy.spatial import ConvexHull  # First thought this warning was caused by a pythonpath problem, but it seems more likely that the warning is caused by scipy import hackery. pylint: disable=no-name-in-module
+    from blmath.geometry import Polyline
+    from blmath.geometry.transform.rotation import rotate_to_xz_plane
+    from blmath.geometry.transform.translation import translation
+
+    if len(polyline.v) <= 1:
+        return polyline
+
+    rotated, R, p0 = rotate_to_xz_plane(polyline.v)
+
+    if flatten:
+        rotated[:, 1] = np.mean(rotated[:, 1])
+
+    projected = rotated[:, [0, 2]]
+
+    hull_vertices = ConvexHull(projected).vertices
+    rotated_hull = rotated[hull_vertices]
+
+    R_inv = np.array(np.mat(R.T).I)
+    inv_rotated = np.dot(rotated_hull, R_inv)
+
+    if p0 is not None:
+        curve_hull = translation(np.array(inv_rotated), -p0)[0]
+    else:
+        curve_hull = inv_rotated
+
+    result = Polyline(v=curve_hull, closed=polyline.closed)
+
+    if want_vertices:
+        return result, hull_vertices
+    else:
+        return result