prepared repo for NASA branch

setup.py

@@ -13,17 +13,26 @@
 from setuptools import setup
 
 desc = 'Microstructure modeling, mesh generation, analysis, and visualization.'
+vs_fname = ('src/microstructpy', '_vs.py')
 
 
 def read(*fname):
     return open(join(dirname(__file__), *fname)).read()
 
 
 def find_version(*fname):
+    ver_str = ''
     for line in read(*fname).split('\n'):
         if line.startswith('__version__') and '=' in line:
-            return line.split('=')[-1].strip().strip('\"').strip('\'')
-    return ''
+            ver_str = line.split('=')[-1].strip().strip('\"').strip('\'')
+            break
+
+    for line in read(*vs_fname).split('\n'):
+        if '+' in line:
+            tag = line.split('+')[-1].strip().strip('\"').strip('\'')
+            ver_str += '+' + tag
+            break
+    return ver_str
 
 
 setup(

src/microstructpy/__init__.py

@@ -3,5 +3,7 @@
 import microstructpy.meshing
 import microstructpy.seeding
 import microstructpy.verification
+from microstructpy._vs import _ver_str
 
 __version__ = '1.1.1'
+__version__ = _ver_str(__version__)

src/microstructpy/_vs.py

@@ -0,0 +1,2 @@
+def _ver_str(version_string):
+    return version_string

src/microstructpy/geometry/_ellipse_best_fit.py

@@ -0,0 +1,68 @@
+import numpy as np
+
+
+def _best_fit(points, ellipse):
+    # Unpack the input points
+    pts = np.array(points, dtype='float')
+    x, y = pts.T
+
+    # Quadratic part of design matrix
+    D1 = np.mat(np.vstack([x * x, x * y, y * y])).T
+
+    # Linear part of design matrix
+    D2 = np.mat(np.vstack([x, y, np.ones(len(x))])).T
+
+    # Scatter matrix
+    S1 = D1.T * D1
+    S2 = D1.T * D2
+    S3 = D2.T * D2
+
+    # Constraint matrix
+    C1inv = np.mat([[0, 0, 0.5], [0, -1, 0], [0.5, 0, 0]])
+
+    # Reduced scatter matrix
+    M = C1inv * (S1 - S2 * S3.I * S2.T)
+
+    # Find eigenvalues
+    _, evec = np.linalg.eig(M)
+
+    # Mask
+    cond = 4 * np.multiply(evec[0, :], evec[2, :])
+    cond -= np.multiply(evec[1, :], evec[1, :])
+    a1 = evec[:, np.nonzero(cond.A > 0)[1]]
+
+    a2 = -S3.I * S2.T * a1
+
+    # Coefficients
+    a = a1[0, 0]
+    b = 0.5 * a1[1, 0]
+    c = a1[2, 0]
+
+    d = 0.5 * a2[0, 0]
+    f = 0.5 * a2[1, 0]
+    g = a2[2, 0]
+
+    # Center of ellipse
+    k = b * b - a * c
+    xc = (c * d - b * f) / k
+    yc = (a * f - b * d) / k
+
+    # Semi-axes lengths
+    numer = a * f * f
+    numer += c * d * d
+    numer += g * b * b
+    numer -= 2 * b * d * f
+    numer -= a * c * g
+    numer *= 2
+
+    tan2 = 2 * b / (a - c)
+    sq_val = np.sqrt(1 + tan2 * tan2)
+    denom1 = k * ((c - a) * sq_val - (c + a))
+    denom2 = k * ((a - c) * sq_val - (c + a))
+    width = np.sqrt(numer / denom1)
+    height = np.sqrt(numer / denom2)
+
+    # Angle of rotation
+    phi = 0.5 * np.arctan(tan2)
+
+    return width, height, phi, xc, yc

src/microstructpy/geometry/ellipse.py

@@ -9,6 +9,8 @@
 import numpy as np
 from matplotlib import patches
 
+from microstructpy.geometry._ellipse_best_fit import _best_fit
+
 __author__ = 'Kenneth (Kip) Hart'
 
 
@@ -47,6 +49,23 @@ class Ellipse(object):
     # Constructor                                                             #
     # ----------------------------------------------------------------------- #
     def __init__(self, **kwargs):
+        # Check Values
+        for kw in ('a', 'b', 'size', 'aspect_ratio'):
+            if kw in kwargs and kwargs[kw] <= 0:
+                raise ValueError(kw + ' should be positive.')
+        if 'axes' in kwargs:
+            for i, ax in kwargs['axes']:
+                if ax <= 0:
+                    raise ValueError('axes[{}] should be positive'.format(i))
+
+        for kw in ['matrix', 'orientation']:
+            if kw in kwargs:
+                m = np.array(kwargs[kw])
+                if m.shape != (2, 2):
+                    raise ValueError(kw + ' should be 2x2.')
+                if not np.all(np.isclose(m * m.T, np.eye(2))):
+                    raise ValueError(kw + ' should be orthonormal.')
+
         # position
         if 'center' in kwargs:
             self.center = kwargs['center']
@@ -173,68 +192,14 @@ def best_fit(self, points):
           (http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.1.7559&rep=rep1&type=pdf)
 
         """  # NOQA: E501
-        # Unpack the input points
-        pts = np.array(points, dtype='float')
-        x, y = pts.T
-
-        # Quadratic part of design matrix
-        D1 = np.mat(np.vstack([x * x, x * y, y * y])).T
-
-        # Linear part of design matrix
-        D2 = np.mat(np.vstack([x, y, np.ones(len(x))])).T
-
-        # Scatter matrix
-        S1 = D1.T * D1
-        S2 = D1.T * D2
-        S3 = D2.T * D2
-
-        # Constraint matrix
-        C1inv = np.mat([[0, 0, 0.5], [0, -1, 0], [0.5, 0, 0]])
-
-        # Reduced scatter matrix
-        M = C1inv * (S1 - S2 * S3.I * S2.T)
-
-        # Find eigenvalues
-        _, evec = np.linalg.eig(M)
-
-        # Mask
-        cond = 4 * np.multiply(evec[0, :], evec[2, :])
-        cond -= np.multiply(evec[1, :], evec[1, :])
-        a1 = evec[:, np.nonzero(cond.A > 0)[1]]
-
-        a2 = -S3.I * S2.T * a1
-
-        # Coefficients
-        a = a1[0, 0]
-        b = 0.5 * a1[1, 0]
-        c = a1[2, 0]
-
-        d = 0.5 * a2[0, 0]
-        f = 0.5 * a2[1, 0]
-        g = a2[2, 0]
-
-        # Center of ellipse
-        k = b * b - a * c
-        xc = (c * d - b * f) / k
-        yc = (a * f - b * d) / k
-
-        # Semi-axes lengths
-        numer = a * f * f
-        numer += c * d * d
-        numer += g * b * b
-        numer -= 2 * b * d * f
-        numer -= a * c * g
-        numer *= 2
-
-        tan2 = 2 * b / (a - c)
-        sq_val = np.sqrt(1 + tan2 * tan2)
-        denom1 = k * ((c - a) * sq_val - (c + a))
-        denom2 = k * ((a - c) * sq_val - (c + a))
-        width = np.sqrt(numer / denom1)
-        height = np.sqrt(numer / denom2)
-
-        # Angle of rotation
-        phi = 0.5 * np.arctan(tan2)
+
+        pts = np.array(points)
+        pt_cen = pts.mean(axis=0)
+        pts -= pt_cen.reshape(1, -1)
+
+        width, height, phi, xc, yc = _best_fit(pts, self)
+        xc += pt_cen[0]
+        yc += pt_cen[1]
 
         # Find pair closest to self
         s = np.sin(phi)

src/microstructpy/seeding/seed.py

@@ -66,7 +66,7 @@ def __init__(self, seed_geometry, phase=0, breakdown=None, position=None):
         self.geometry = seed_geometry
         self.phase = phase
 
-        if position is None:
+        if position is None and self.geometry is not None:
             self.position = [0 for _ in range(self.geometry.n_dim)]
         else:
             self.position = position
@@ -121,15 +121,22 @@ def factory(cls, seed_type, phase=0, breakdown=None, position=None,
         assert type(seed_type) is str
         seed_type = seed_type.strip().lower()
 
-        n_dim = geometry.factory(seed_type).n_dim
+        if seed_type == 'nonetype':
+            n_dim = 0
+        else:
+            n_dim = geometry.factory(seed_type).n_dim
         if 'volume' in kwargs:
             if n_dim == 2:
                 size = 2 * np.sqrt(kwargs['volume'] / np.pi)
             else:
                 size = 2 * np.cbrt(3 * kwargs['volume'] / (4 * np.pi))
             kwargs['size'] = size
 
-        geom = geometry.factory(seed_type, **kwargs)
+        # Catch NoneType geometries
+        if seed_type == 'nonetype':
+            geom = None
+        else:
+            geom = geometry.factory(seed_type, **kwargs)
 
         if breakdown is None:
             if seed_type in ('circle', 'sphere'):
@@ -162,10 +169,14 @@ def from_str(cls, seed_str):
         # Convert to dictionary
         str_dict = {}
         for line in seed_str.strip().split('\n'):
-            k_str, v_str = line.split(':')
-            k = k_str.strip().lower().replace(' ', '_')
-            v = _misc.from_str(v_str)
-            str_dict[k] = v
+            try:
+                k_str, v_str = line.split(':')
+            except ValueError:
+                continue
+            else:
+                k = k_str.strip().lower().replace(' ', '_')
+                v = _misc.from_str(v_str)
+                str_dict[k] = v
 
         # Extract seed type, phase, and breakdown
         seed_type = str_dict['geometry']
@@ -220,6 +231,11 @@ def __repr__(self):
     # Comparison Functions                                                    #
     # ----------------------------------------------------------------------- #
     def __lt__(self, seed):
+        if self.geometry is None:
+            return True
+        if seed.geometry is None:
+            return False
+
         a_str = 'Seeds are not the same dimension.'
         assert self.geometry.n_dim == seed.geometry.n_dim, a_str
         if self.geometry.n_dim == 2:
@@ -300,6 +316,8 @@ def position(self, pos):
     @property
     def volume(self):
         """float: The area (2D) or volume (3D) of the seed"""
+        if self.geometry is None:
+            return 0
         if self.geometry.n_dim == 2:
             return self.geometry.area
         else:
@@ -311,6 +329,8 @@ def volume(self):
     @property
     def limits(self):
         """list: The (lower, upper) bounds of the seed"""
+        if self.geometry is None:
+            return []
         return self.geometry.limits
 
     # ----------------------------------------------------------------------- #
@@ -328,7 +348,8 @@ def plot(self, **kwargs):
             **kwargs: Plotting keyword arguments.
 
         """
-        self.geometry.plot(**kwargs)
+        if self.geometry is not None:
+            self.geometry.plot(**kwargs)
 
     # ----------------------------------------------------------------------- #
     # Plot Breakdown                                                          #
@@ -345,7 +366,7 @@ def plot_breakdown(self, **kwargs):
 
         """
 
-        n = self.geometry.n_dim
+        n = len(self.breakdown[0]) - 1
         if n == 2:
             pc = [patches.Circle([x, y], r) for x, y, r in self.breakdown]
             coll = collections.PatchCollection(pc, **kwargs)

src/microstructpy/seeding/seedlist.py

@@ -454,8 +454,11 @@ def plot(self, **kwargs):
                         val_list.append(val)
                         pc_kwargs[key] = val_list
 
+                elif geom_name == 'nonetype':
+                    pass
+
                 else:
-                    e_str = 'Cannot plot groups of ' + str(seed.seed_type)
+                    e_str = 'Cannot plot groups of ' + geom_name
                     e_str += ' yet.'
                     raise NotImplementedError(e_str)