modifications to findorientations + bug fixes

Author: Joel Bernier

hexrd/config/findorientations.py

@@ -70,7 +70,7 @@ class ClusteringConfig(Config):
     def algorithm(self):
         key = 'find_orientations:clustering:algorithm'
         choices = ['dbscan', 'ort-dbscan', 'sph-dbscan', 'fclusterdata']
-        temp = self._cfg.get(key, 'ort-dbscan').lower()
+        temp = self._cfg.get(key, 'dbscan').lower()
         if temp in choices:
             return temp
         raise RuntimeError(

hexrd/findorientations.py

@@ -173,10 +173,10 @@ def quat_distance(x, y):
         qfib_r = qfib[:, np.array(compl) > min_compl]
 
         num_ors = qfib_r.shape[1]
-        
+
         if num_ors > 25000:
             if algorithm == 'sph-dbscan' or algorithm == 'fclusterdata':
-                logger.info("defaulting to orthographic DBSCAN")
+                logger.info("falling back to euclidean DBSCAN")
                 algorithm = 'ort-dbscan'
             #raise RuntimeError, \
             #    "Requested clustering of %d orientations, which would be too slow!" %qfib_r.shape[1]
@@ -191,18 +191,19 @@ def quat_distance(x, y):
             logger.warning(
                 "sklearn >= 0.14 required for dbscan; using fclusterdata"
                 )
-                
+
         if algorithm == 'dbscan' or algorithm == 'ort-dbscan' or algorithm == 'sph-dbscan':
             # munge min_samples according to options
             if min_samples is None or cfg.find_orientations.use_quaternion_grid is not None:
                 min_samples = 1
-            
+
             if algorithm == 'sph-dbscan':
+                logger.info("using spherical DBSCAN")
                 # compute distance matrix
                 pdist = pairwise_distances(
                     qfib_r.T, metric=quat_distance, n_jobs=1
                     )
-    
+
                 # run dbscan
                 core_samples, labels = dbscan(
                     pdist,
@@ -212,11 +213,13 @@ def quat_distance(x, y):
                     )
             else:
                 if algorithm == 'ort-dbscan':
+                    logger.info("using euclidean orthographic DBSCAN")
                     pts = qfib_r[1:, :].T
-                    eps = 0.9*np.radians(cl_radius)
+                    eps = 0.25*np.radians(cl_radius)
                 else:
+                    logger.info("using euclidean DBSCAN")
                     pts = qfib_r.T
-                    eps = np.radians(cl_radius)
+                    eps = 0.5*np.radians(cl_radius)
 
                 # run dbscan
                 core_samples, labels = dbscan(
@@ -226,19 +229,20 @@ def quat_distance(x, y):
                     metric='minkowski', p=2,
                     )
 
-            # extract cluster labels    
+            # extract cluster labels
             cl = np.array(labels, dtype=int) # convert to array
             noise_points = cl == -1 # index for marking noise
             cl += 1 # move index to 1-based instead of 0
             cl[noise_points] = -1 # re-mark noise as -1
-            logger.info("dbscan found %d noise points", sum(noise_points))     
+            logger.info("dbscan found %d noise points", sum(noise_points))
         elif algorithm == 'fclusterdata':
+            logger.info("using spherical fclusetrdata")
             cl = cluster.hierarchy.fclusterdata(
                 qfib_r.T,
                 np.radians(cl_radius),
                 criterion='distance',
                 metric=quat_distance
-                )      
+                )
         else:
             raise RuntimeError(
                 "Clustering algorithm %s not recognized" % algorithm
@@ -249,26 +253,18 @@ def quat_distance(x, y):
             nblobs = len(np.unique(cl)) - 1
         else:
             nblobs = len(np.unique(cl))
-        
+
         #import pdb; pdb.set_trace()
-        
+
         """ PERFORM AVERAGING TO GET CLUSTER CENTROIDS """
         qbar = np.zeros((4, nblobs))
-        if algorithm == 'sph-dbscan' or algorithm == 'fclusterdata':
-            # here clusters can be split across fr
-            for i in range(nblobs):
-                npts = sum(cl == i + 1) 
-                qbar[:, i] = rot.quatAverageCluster(
-                    qfib_r[:, cl == i + 1].reshape(4, npts), qsym
-                    ).flatten()
-                pass
-        else:
-            # here clusters are ompact by construction
-            for i in range(nblobs):
-                qbar[:, i] = np.average(np.atleast_2d(qfib_r[:, cl == i + 1]), axis=1)
-                pass
-            qbar = sym.toFundamentalRegion(mutil.unitVector(qbar), crysSym=qsym)
-
+        for i in range(nblobs):
+            npts = sum(cl == i + 1)
+            qbar[:, i] = rot.quatAverageCluster(
+                qfib_r[:, cl == i + 1].reshape(4, npts), qsym
+            ).flatten()
+            pass
+        pass
     logger.info("clustering took %f seconds", time.clock() - start)
     logger.info(
         "Found %d orientation clusters with >=%.1f%% completeness"
@@ -377,6 +373,8 @@ def find_orientations(cfg, hkls=None, clean=False, profile=False):
     else:
         distortion = None
 
+    min_compl = cfg.find_orientations.clustering.completeness
+
     # start logger
     logger.info("beginning analysis '%s'", cfg.analysis_name)
 
@@ -477,8 +475,9 @@ def find_orientations(cfg, hkls=None, clean=False, profile=False):
             refl_per_grain[i] = len(sim_results[0])
             num_seed_refls[i] = np.sum([sum(sim_results[0] == hkl_id) for hkl_id in hkl_ids])
             pass
+        #min_samples = 2
         min_samples = max(
-            np.floor(cfg.find_orientations.clustering.completeness*np.average(num_seed_refls)),
+            int(np.floor(0.5*min_compl*min(num_seed_refls))),
             2
             )
         mean_rpg = int(np.round(np.average(refl_per_grain)))

hexrd/xrd/crystallography.py

@@ -863,15 +863,15 @@ def getMergedRanges(self, cullDupl=False):
 
         # if you end exlcusions in a doublet (or multiple close rings)
         # then this will 'fail'.  May need to revisit...
-        nonoverlapNexts = np.hstack((tThRanges[:-1,1] < tThRanges[1:,0], True))
+        nonoverlapNexts = num.hstack((tThRanges[:-1,1] < tThRanges[1:,0], True))
         iHKLLists = []
         mergedRanges = []
         hklsCur = []
         tThLoIdx = 0
-        tThLoCur = tThHiCur = 0.
+        tThHiCur = 0.
         for iHKL, nonoverlapNext in enumerate(nonoverlapNexts):
             print tThLoIdx
-            ThHi = tThRanges[iHKL, -1]
+            tThHi = tThRanges[iHKL, -1]
             if not nonoverlapNext:
                 if cullDupl and abs(tThs[iHKL] - tThs[iHKL+1]) < sqrt_epsf:
                   continue

hexrd/xrd/rotations.py

@@ -352,18 +352,10 @@ def quatAverageCluster(q_in, qsym):
             q_bar = quatProduct(q_in[:, 0].reshape(4, 1),
                                 quatOfExpMap(0.5*ma*unitVector(mq[1:].reshape(3, 1))))
     else:
-        # use first quat as initial guess
-        phi = 2. * arccos(q_in[0, 0])
-        if phi <= finfo(float).eps:
-            x0 = zeros(3)
-        else:
-            n = unitVector(q_in[1:, 0].reshape(3, 1))
-            x0 = phi*n.flatten()
-
         # first drag to origin using first quat (arb!)
         q0 = q_in[:, 0].reshape(4, 1)
         qrot = dot(
-            quatProductMatrix( invertQuat( q0 ), mult='right' ),
+            quatProductMatrix( invertQuat( q0 ), mult='left' ),
             q_in)
 
         # second, re-cast to FR
@@ -374,7 +366,7 @@ def quatAverageCluster(q_in, qsym):
 
         # unrotate!
         q_bar = dot(
-            quatProductMatrix( q0, mult='right' ),
+            quatProductMatrix( q0, mult='left' ),
             q_bar)
 
         # re-map

hexrd/xrd/xrdutil.py

@@ -3170,53 +3170,6 @@ def validateAngleRanges(angList, startAngs, stopAngs, ccw=True):
                     reflInRange = reflInRange | num.logical_and(zStart >= 0, zStop <= 0)
     return reflInRange
 
-def tVec_d_from_old_parfile(old_par, detOrigin):
-    beamXYD = old_par[:3, 0]
-    rMat_d  = xf.makeDetectorRotMat(old_par[3:6, 0])
-    bVec_ref = num.c_[0., 0., -1.].T
-    args=(rMat_d, beamXYD, detOrigin, bVec_ref)
-    tvd_xy = opt.leastsq(objFun_tVec_d, -beamXYD[:2], args=args)[0]
-    return num.hstack([tvd_xy, -beamXYD[2]]).reshape(3, 1)
-
-def objFun_tVec_d(tvd_xy, rMat_d, beamXYD, detOrigin, bHat_l):
-    """
-    """
-    xformed_xy = beamXYD[:2] - detOrigin
-    tVec_d = num.hstack([tvd_xy, -beamXYD[2]]).T
-    n_d    = rMat_d[:, 2]
-
-    bVec_l = (num.dot(n_d, tVec_d) / num.dot(n_d, bHat_l)) * bHat_l
-    bVec_d = num.hstack([xformed_xy, 0.]).T
-
-    return num.dot(rMat_d, bVec_d).flatten() + tVec_d.flatten() - bVec_l.flatten()
-
-def beamXYD_from_tVec_d(rMat_d, tVec_d, bVec_ref, detOrigin):
-    # calculate beam position
-    Zd_l = num.dot(rMat_d, num.c_[0, 0, 1].T)
-    bScl = num.dot(Zd_l.T, tVec_d) / num.dot(Zd_l.T, bVec_ref)
-    beamPos_l = bScl*bVec_ref
-    return num.dot(rMat_d.T, beamPos_l - tVec_d) + num.hstack([detOrigin, -tVec_d[2]]).reshape(3, 1)
-
-def write_old_parfile(filename, results):
-    if isinstance(filename, file):
-        fid = filename
-    elif isinstance(filename, str) or isinstance(filename, unicode):
-        fid = open(filename, 'w')
-        pass
-    rMat_d = xf.makeDetectorRotMat(results['tiltAngles'])
-    tVec_d = results['tVec_d'] - results['tVec_s']
-    beamXYD = beamXYD_from_tVec_d(rMat_d, tVec_d, bVec_ref, detOrigin)
-    det_plist = num.zeros(12)
-    det_plist[:3]  = beamXYD.flatten()
-    det_plist[3:6] = results['tiltAngles']
-    det_plist[6:]  = results['dParams']
-    print >> fid, "# DETECTOR PARAMETERS (from new geometry model fit)"
-    print >> fid, "# \n# <class 'hexrd.xrd.detector.DetectorGeomGE'>\n#"
-    for i in range(len(det_plist)):
-        print >> fid, "%1.8e\t%d" % (det_plist[i], 0)
-    fid.close()
-    return
-
 def simulateOmeEtaMaps(omeEdges, etaEdges, planeData, expMaps,
                        chi=0.,
                        etaTol=None, omeTol=None,
@@ -3429,7 +3382,7 @@ def simulateGVecs(pd, detector_params, grain_params,
                   pixel_pitch=(0.2, 0.2),
                   distortion=(dFunc_ref, dParams_ref)):
     """
-    returns valid_hkl, valid_ang, valid_xy, ang_ps
+    returns valid_ids, valid_hkl, valid_ang, valid_xy, ang_ps
 
     panel_dims are [(xmin, ymin), (xmax, ymax)] in mm
 
@@ -3506,17 +3459,19 @@ def simulateGVecs(pd, detector_params, grain_params,
             det_xy[:, 1] >= panel_dims[0][1], 
             det_xy[:, 1] <= panel_dims[1][1]
             )
-        on_panel   = num.logical_and(on_panel_x, on_panel_y)
+        on_panel = num.logical_and(on_panel_x, on_panel_y)
         #
-        valid_ang = allAngs[on_panel, :]
+        op_idx = num.where(on_panel)[0]
+        #
+        valid_ang = allAngs[op_idx, :]
         valid_ang[:, 2] = xf.mapAngle(valid_ang[:, 2], ome_period)
-        valid_ids = allHKLs[on_panel, 0]
-        valid_hkl = allHKLs[on_panel, 1:]
-        valid_xy  = det_xy[on_panel, :]
-        ang_ps    = angularPixelSize(valid_xy, pixel_pitch,
-                                     rMat_d, rMat_s,
-                                     tVec_d, tVec_s, tVec_c,
-                                     distortion=distortion)
+        valid_ids = allHKLs[op_idx, 0]
+        valid_hkl = allHKLs[op_idx, 1:]
+        valid_xy = det_xy[op_idx, :]
+        ang_ps = angularPixelSize(valid_xy, pixel_pitch,
+                                  rMat_d, rMat_s,
+                                  tVec_d, tVec_s, tVec_c,
+                                  distortion=distortion)
 
     return valid_ids, valid_hkl, valid_ang, valid_xy, ang_ps
 
@@ -3575,7 +3530,6 @@ def simulateLauePattern(hkls, bMat,
     LOOP OVER GRAINS
     """   
 
-    
     for iG, gp in enumerate(grain_params):
         rmat_c = xfcapi.makeRotMatOfExpMap(gp[:3])
         tvec_c = gp[3:6].reshape(3, 1)