Minor simplification and improvement to the final check for diverging/invalid coordinates compared to r31125. Improvements/corrections in code formatting, comments, and help.

git-svn-id: http://svn.stsci.edu/svn/ssb/stsci_python/stwcs/trunk@31126 fe389314-cf27-0410-b35b-8c050e845b92

1 files changed, 50 insertions, 43 deletions

diff --git a/lib/stwcs/wcsutil/hstwcs.py b/lib/stwcs/wcsutil/hstwcs.py
index 8f66283..f48741d 100644
--- a/lib/stwcs/wcsutil/hstwcs.py
+++ b/lib/stwcs/wcsutil/hstwcs.py
@@ -446,7 +446,8 @@ class HSTWCS(WCS):
 
     def all_sky2pix(self, *args, **kwargs):
         """
-        all_sky2pix(*arg, accuracy=1.0e-4, maxiter=20, adaptive=False, quiet=False)
+        all_sky2pix(*arg, accuracy=1.0e-4, maxiter=20, adaptive=False, \
+detect_divergence=True, quiet=False)
 
         Performs full inverse transformation using iterative solution
         on full forward transformation with complete distortion model.
@@ -456,7 +457,7 @@ class HSTWCS(WCS):
         accuracy : float, optional (Default = 1.0e-4)
             Required accuracy of the solution. Iteration terminates when the
             correction to the solution found during the previous iteration
-            is smaller (in the sence of the L2 norm) than `tolerance`\ .
+            is smaller (in the sence of the L2 norm) than `accuracy`\ .
 
         maxiter : int, optional (Default = 20)
             Maximum number of iterations allowed to reach the solution.
@@ -621,7 +622,7 @@ accuracy after 3 iterations.
         [[  5.52645241 -72.05171776]
          [  7.15979392 -70.81405561]
          [  5.52653313 -72.05170814]]
-         
+
         >>> try:
         >>>   xy = w.all_sky2pix(divradec,1, maxiter=20, accuracy=1.0e-4, \
 adaptive=False, detect_divergence=True, quiet=False)
@@ -647,7 +648,7 @@ adaptive=False, detect_divergence=True, quiet=False)
           File "<console>", line 8, in <module>
         NoConvergence: 'HSTWCS.all_sky2pix' failed to converge to the requested accuracy.
         After 5 iterations, the solution is diverging at least for one input point.
-        
+
         >>> try:
         >>>   xy = w.all_sky2pix(divradec,1, maxiter=20, accuracy=1.0e-4, \
 adaptive=False, detect_divergence=False, quiet=False)
@@ -675,6 +676,9 @@ adaptive=False, detect_divergence=False, quiet=False)
         After 20 iterations, the solution is diverging at least for one input point.
 
         """
+        #####################################################################
+        ##                     PROCESS ARGUMENTS:                          ##
+        #####################################################################
         nargs = len(args)
 
         if nargs == 3:
@@ -710,13 +714,15 @@ adaptive=False, detect_divergence=False, quiet=False)
         detect_divergence = kwargs.pop('detect_divergence', True)
         quiet             = kwargs.pop('quiet', False)
 
-        # initialize iterative process:
+        #####################################################################
+        ##                INITIALIZE ITERATIVE PROCESS:                    ##
+        #####################################################################
         x0, y0 = self.wcs_sky2pix(ra, dec, origin) # <-- initial approximation
                                                   #     (WCS based only)
 
-        # see if iterative solution is required (when any of the
+        # see if an iterative solution is required (when any of the
         # non-CD-matrix corrections are present). If not required
-        # return initial approximation (xy0).
+        # return initial approximation (x0, y0).
         if self.sip is None and \
            self.cpdis1 is None and self.cpdis2 is None and \
            self.det2im1 is None and self.det2im2 is None:
@@ -747,19 +753,22 @@ adaptive=False, detect_divergence=False, quiet=False)
         dn2prev   = dx**2+dy**2
         dn2       = dn2prev
 
-        # process all coordinates simultaneously:
+        # prepare for iterative process
         iterlist  = range(1, maxiter+1)
         accuracy2 = accuracy**2
         ind       = None
         inddiv    = None
 
         npts      = x.shape[0]
-        
-        # turn off numpy runtime warnings for 'invalid' and 'over':        
+
+        # turn off numpy runtime warnings for 'invalid' and 'over':
         old_invalid = np.geterr()['invalid']
         old_over    = np.geterr()['over']
         np.seterr(invalid = 'ignore', over = 'ignore')
 
+        #####################################################################
+        ##                     NON-ADAPTIVE ITERATIONS:                    ##
+        #####################################################################
         if not adaptive:
             for k in iterlist:
                 # check convergence:
@@ -803,8 +812,10 @@ adaptive=False, detect_divergence=False, quiet=False)
                 x -= dx
                 y -= dy
 
+        #####################################################################
+        ##                      ADAPTIVE ITERATIONS:                       ##
+        #####################################################################
         if adaptive:
-
             if ind is None:
                 ind = np.asarray(range(npts), dtype=np.int64)
 
@@ -828,50 +839,43 @@ adaptive=False, detect_divergence=False, quiet=False)
                 # update indices of elements that still need correction:
                 if detect_divergence:
                     ind, = np.where((dn2 >= accuracy2) & (dn2 <= dn2prev))
+                    #ind = ind[np.where((dn2[ind] >= accuracy2) & (dn2[ind] <= dn2prev))]
                     dn2prev[ind] = dn2[ind]
                 else:
                     ind, = np.where(dn2 >= accuracy2)
-                    #ind = ind[np.where(dx[ind]**2+dy[ind]**2 >= accuracy2)]
+                    #ind = ind[np.where(dn2[ind] >= accuracy2)]
 
                 # apply correction:
                 x[ind] -= dx[ind]
                 y[ind] -= dy[ind]
 
-        #FINAL DETECTION OF INVALID, DIVERGING, AND FAILED-TO-CONVERGE POINTS:
+        #####################################################################
+        ##         FINAL DETECTION OF INVALID, DIVERGING,                  ##
+        ##         AND FAILED-TO-CONVERGE POINTS                           ##
+        #####################################################################
+        # Identify diverging and/or invalid points:
         invalid = (((~np.isfinite(y)) | (~np.isfinite(x)) | \
                     (~np.isfinite(dn2))) & \
                    (np.isfinite(ra)) &  (np.isfinite(dec)))
-        if detect_divergence:
-            inddiv, = np.where(((dn2 >= accuracy2) & (dn2 > dn2prev)) |   \
-                               invalid)
-            # identify diverging and/or invalid points:
-            if inddiv.shape[0] == 0:
-                inddiv = None
-            # identify points that did not converge within
-            # 'maxiter' iterations:
-            if k >= maxiter:
-                ind, = np.where((dn2 >= accuracy2) & (dn2 <= dn2prev) & \
-                                (~invalid))
-                if ind.shape[0] == 0:
-                    ind = None
-            else:
+        # When detect_divergence==False, dn2prev is outdated (it is the
+        # norm^2 of the very first correction). Still better than nothing...
+        inddiv, = np.where(((dn2 >= accuracy2) & (dn2 > dn2prev)) | invalid)
+        if inddiv.shape[0] == 0:
+            inddiv = None
+        # identify points that did not converge within
+        # 'maxiter' iterations:
+        if k >= maxiter:
+            ind,= np.where((dn2 >= accuracy2) & (dn2 <= dn2prev) & (~invalid))
+            if ind.shape[0] == 0:
                 ind = None
-
         else:
-            # dn2prev is not available in this case.
-            inddiv, = np.where(invalid)
-            # identify diverging and/or invalid points:
-            if inddiv.shape[0] == 0:
-                inddiv = None
-            # identify points that did not converge within
-            # 'maxiter' iterations:
-            if k >= maxiter:
-                ind, = np.where((dn2 >= accuracy2) & (~invalid))
-                if ind.shape[0] == 0:
-                    ind = None
-            else:
-                ind = None
+            ind = None
 
+        #####################################################################
+        ##      RAISE EXCEPTION IF DIVERGING OR TOO SLOWLY CONVERGING      ##
+        ##      DATA POINTS HAVE BEEN DETECTED:                            ##
+        #####################################################################
+        # raise exception if diverging or too slowly converging
         if (ind is not None or inddiv is not None) and not quiet:
             if vect1D:
                 sol  = [x, y]
@@ -882,8 +886,8 @@ adaptive=False, detect_divergence=False, quiet=False)
 
             # restore previous numpy error settings:
             np.seterr(invalid = old_invalid, over = old_over)
-            
-            if inddiv is None:                
+
+            if inddiv is None:
                 raise NoConvergence("'HSTWCS.all_sky2pix' failed to "        \
                     "converge to the requested accuracy after {:d} "         \
                     "iterations.".format(k), best_solution = sol,            \
@@ -898,6 +902,9 @@ adaptive=False, detect_divergence=False, quiet=False)
                     accuracy = err, niter = k, failed2converge = ind,        \
                     divergent = inddiv)
 
+        #####################################################################
+        ##             FINALIZE AND FORMAT DATA FOR RETURN:                ##
+        #####################################################################
         # restore previous numpy error settings:
         np.seterr(invalid = old_invalid, over = old_over)