This revision improves upon changes made to the `all_sky2pix` in r31122 and r31123 that address issues described in ticket #1121 and ticket #1122, by implementing a more robust handling of the divergent solutions. The new parameter `detect_divergence` will allow to detect divergent solutions at each iteration. Now, even if `adaptive` is set to `False`, assuming `detect_divergence` is `True`, when divergent solutions are detected `all_sky2pix` will automatically switch from the highly efficient fully vectorized algorithm to a less efficient (~30%) (still vectorized) but adaptive algorithm.

Updated comments and examples to `all_sky2pix` to reflect the changes.



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

1 files changed, 271 insertions, 83 deletions

diff --git a/lib/stwcs/wcsutil/hstwcs.py b/lib/stwcs/wcsutil/hstwcs.py
index e6c812d..8f66283 100644
--- a/lib/stwcs/wcsutil/hstwcs.py
+++ b/lib/stwcs/wcsutil/hstwcs.py
@@ -1,6 +1,6 @@
 from __future__ import division # confidence high
 
-import os.path
+import os
 from pywcs import WCS
 import pyfits
 import instruments
@@ -79,12 +79,12 @@ class NoConvergence(Exception):
         solution appears to be divergent. If the solution does not diverge,
         `divergent` will be set to `None`.
 
-    nonconvergent : None, numpy.array
+    failed2converge : None, numpy.array
         Indices of the points in :py:attr:`best_solution` array for which the
         solution failed to converge within the specified maximum number
         of iterations. If there are no non-converging poits (i.e., if
         the required accuracy has been achieved for all points) then
-        `nonconvergent` will be set to `None`.
+        `failed2converge` will be set to `None`.
 
     """
     def __init__(self, *args, **kwargs):
@@ -94,7 +94,7 @@ class NoConvergence(Exception):
         self.accuracy       = kwargs.pop('accuracy', None)
         self.niter          = kwargs.pop('niter', None)
         self.divergent      = kwargs.pop('divergent', None)
-        self.nonconvergent  = kwargs.pop('nonconvergent', None)
+        self.failed2converge= kwargs.pop('failed2converge', None)
 
 
 #
@@ -454,7 +454,9 @@ class HSTWCS(WCS):
         Parameters
         ----------
         accuracy : float, optional (Default = 1.0e-4)
-            Required accuracy of the solution.
+            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`\ .
 
         maxiter : int, optional (Default = 20)
             Maximum number of iterations allowed to reach the solution.
@@ -462,25 +464,94 @@ class HSTWCS(WCS):
         adaptive : bool, optional (Default = False)
             Specifies whether to adaptively select only points that did not
             converge to a solution whithin the required accuracy for the
-            next iteration. Default is recommended for HST instruments.
+            next iteration. Default is recommended for HST as well as most
+            other instruments.
 
             .. note::
-               The :py:meth:`all_sky2pix` uses a vectorized implementation of
-               the method of consecutive approximations (see `Notes` section
-               below) in which it iterates over *all* input poits *regardless*
-               untill the required accuracy has been reached for *all* input
-               points. In some cases it may be possible that *almost all* points
-               have reached the required accuracy but there are only a few
-               of input data points for which additional iterations may be
-               needed. In this situation it may be advantageous to set
-               `adaptive` = `True`\ in which case :py:meth:`all_sky2pix` will
-               continue iterating *only* over the points that have not yet
-               converged to the required accuracy. However, for the HST's
-               ACS/WFC detector, which has the strongest distortions of all
-               HST instruments, testing has shown that enabling this option
-               would lead to a 10-30\% penalty in computational time.
-               Therefore, for HST instruments, it is recommended to set
-               `adaptive` = `False`\ .
+               The :py:meth:`all_sky2pix` uses a vectorized implementation
+               of the method of consecutive approximations (see `Notes`
+               section below) in which it iterates over *all* input poits
+               *regardless* until the required accuracy has been reached for
+               *all* input points. In some cases it may be possible that
+               *almost all* points have reached the required accuracy but
+               there are only a few of input data points left for which
+               additional iterations may be needed (this depends mostly on the
+               characteristics of the geometric distortions for a given
+               instrument). In this situation it may be
+               advantageous to set `adaptive` = `True`\ in which
+               case :py:meth:`all_sky2pix` will continue iterating *only* over
+               the points that have not yet converged to the required
+               accuracy. However, for the HST's ACS/WFC detector, which has
+               the strongest distortions of all HST instruments, testing has
+               shown that enabling this option would lead to a about 10-30\%
+               penalty in computational time (depending on specifics of the
+               image, geometric distortions, and number of input points to be
+               converted). Therefore, for HST instruments,
+               it is recommended to set `adaptive` = `False`\ . The only
+               danger in getting this setting wrong will be a performance
+               penalty.
+
+            .. note::
+               When `detect_divergence` is `True`\ , :py:meth:`all_sky2pix` \
+               will automatically switch to the adaptive algorithm once
+               divergence has been detected.
+
+        detect_divergence : bool, optional (Default = True)
+            Specifies whether to perform a more detailed analysis of the
+            convergence to a solution. Normally :py:meth:`all_sky2pix`
+            may not achieve the required accuracy
+            if either the `tolerance` or `maxiter` arguments are too low.
+            However, it may happen that for some geometric distortions
+            the conditions of convergence for the the method of consecutive
+            approximations used by :py:meth:`all_sky2pix` may not be
+            satisfied, in which case consecutive approximations to the
+            solution will diverge regardless of the `tolerance` or `maxiter`
+            settings.
+
+            When `detect_divergence` is `False`\ , these divergent points
+            will be detected as not having achieved the required accuracy
+            (without further details). In addition, if `adaptive` is `False`
+            then the algorithm will not know that the solution (for specific
+            points) is diverging and will continue iterating and trying to
+            "improve" diverging solutions. This may result in NaN or Inf
+            values in the return results (in addition to a performance
+            penalties). Even when `detect_divergence` is
+            `False`\ , :py:meth:`all_sky2pix`\ , at the end of the iterative
+            process, will identify invalid results (NaN or Inf) as "diverging"
+            solutions and will raise :py:class:`NoConvergence` unless
+            the `quiet` parameter is set to `True`\ .
+
+            When `detect_divergence` is `True`\ , :py:meth:`all_sky2pix` will
+            detect points for
+            which current correction to the coordinates is larger than
+            the correction applied during the previous iteration **if** the
+            requested accuracy **has not yet been achieved**\ . In this case,
+            if `adaptive` is `True`, these points will be excluded from
+            further iterations and if `adaptive`
+            is `False`\ , :py:meth:`all_sky2pix` will automatically
+            switch to the adaptive algorithm.
+
+            .. note::
+               When accuracy has been achieved, small increases in
+               current corrections may be possible due to rounding errors
+               (when `adaptive` is `False`\ ) and such increases
+               will be ignored.
+
+            .. note::
+               Setting `detect_divergence` to `True` will incurr about 5-10\%
+               performance penalty (in our testing on ACS/WFC images).
+               Because the benefits of enabling this feature outweigh
+               the small performance penalty, it is recommended to set
+               `detect_divergence` to `True`\ , unless extensive testing
+               of the distortion models for images from specific
+               instruments show a good stability of the numerical method
+               for a wide range of coordinates (even outside the image
+               itself).
+
+            .. note::
+               Indices of the diverging inverse solutions will be reported
+               in the `divergent` attribute of the
+               raised :py:class:`NoConvergence` object.
 
         quiet : bool, optional (Default = False)
             Do not throw :py:class:`NoConvergence` exceptions when the method
@@ -497,9 +568,9 @@ class HSTWCS(WCS):
 
         Notes
         -----
-        Inputs can either be (RA, Dec, origin) or (RADec, origin) where RA and Dec
-        are 1-D arrays/lists of coordinates and RADec is an array/list of pairs
-        of coordinates.
+        Inputs can either be (RA, Dec, origin) or (RADec, origin) where RA
+        and Dec are 1-D arrays/lists of coordinates and RADec is an
+        array/list of pairs of coordinates.
 
         Using the method of consecutive approximations we iterate starting
         with the initial approximation, which is computed using the
@@ -510,8 +581,8 @@ class HSTWCS(WCS):
         highly efficient (>30x) when *all* data points that need to be
         converted from sky coordinates to image coordinates are passed at
         *once*\ . Therefore, it is advisable, whenever possible, to pass
-        as input a long array of all points that need to be converted to
-        :py:meth:`all_sky2pix` instead of calling :py:meth:`all_sky2pix`
+        as input a long array of all points that need to be converted
+        to :py:meth:`all_sky2pix` instead of calling :py:meth:`all_sky2pix`
         for each data point. Also see the note to the `adaptive` parameter.
 
         Examples
@@ -520,6 +591,7 @@ class HSTWCS(WCS):
         >>> hdulist = pyfits.open('j94f05bgq_flt.fits')
         >>> w = stwcs.wcsutil.HSTWCS(hdulist, ext=('sci',1))
         >>> hdulist.close()
+
         >>> ra, dec = w.all_pix2sky([1,2,3],[1,1,1],1); print(ra); print(dec)
         [ 5.52645241  5.52649277  5.52653313]
         [-72.05171776 -72.05171295 -72.05170814]
@@ -537,10 +609,71 @@ class HSTWCS(WCS):
         [[ 1.00000233  0.99999997]
          [ 2.00000232  0.99999997]
          [ 3.00000233  0.99999998]]
-        >>> xy = w.all_sky2pix(radec,1, maxiter=3, accuracy=1.0e-10, quiet=False)
+        >>> xy = w.all_sky2pix(radec,1, maxiter=3, accuracy=1.0e-10, \
+quiet=False)
         NoConvergence: 'HSTWCS.all_sky2pix' failed to converge to requested \
 accuracy after 3 iterations.
 
+        >>>
+        Now try to use some diverging data:
+        >>> divradec = w.all_pix2sky([[1.0,1.0],[10000.0,50000.0],\
+[3.0,1.0]],1); print(divradec)
+        [[  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)
+        >>> except stwcs.wcsutil.hstwcs.NoConvergence as e:
+        >>>   print("Indices of diverging points: {}".format(e.divergent))
+        >>>   print("Indices of poorly converging points: {}".format(e.failed2converge))
+        >>>   print("Best solution: {}".format(e.best_solution))
+        >>>   print("Achieved accuracy: {}".format(e.accuracy))
+        >>>   raise e
+        Indices of diverging points:
+        [1]
+        Indices of poorly converging points:
+        None
+        Best solution:
+        [[  1.00006219e+00   9.99999288e-01]
+         [ -1.99440907e+06   1.44308548e+06]
+         [  3.00006257e+00   9.99999316e-01]]
+        Achieved accuracy:
+        [[  5.98554253e-05   6.79918148e-07]
+         [  8.59514088e+11   6.61703754e+11]
+         [  6.02334592e-05   6.59713067e-07]]
+        Traceback (innermost last):
+          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)
+        >>> except stwcs.wcsutil.hstwcs.NoConvergence as e:
+        >>>   print("Indices of diverging points: {}".format(e.divergent))
+        >>>   print("Indices of poorly converging points: {}".format(e.failed2converge))
+        >>>   print("Best solution: {}".format(e.best_solution))
+        >>>   print("Achieved accuracy: {}".format(e.accuracy))
+        >>>   raise e
+        Indices of diverging points:
+        [1]
+        Indices of poorly converging points:
+        None
+        Best solution:
+        [[  1.   1.]
+         [ nan  nan]
+         [  3.   1.]]
+        Achieved accuracy:
+        [[  0.   0.]
+         [ nan  nan]
+         [  0.   0.]]
+        Traceback (innermost last):
+          File "<console>", line 8, in <module>
+        NoConvergence: 'HSTWCS.all_sky2pix' failed to converge to the requested accuracy.
+        After 20 iterations, the solution is diverging at least for one input point.
+
         """
         nargs = len(args)
 
@@ -571,13 +704,15 @@ accuracy after 3 iterations.
                             .format(nargs))
 
         # process optional arguments:
-        accuracy = kwargs.pop('accuracy', 1.0e-4)
-        maxiter  = kwargs.pop('maxiter', 20)
-        quiet    = kwargs.pop('quiet', False)
-        adaptive = kwargs.pop('adaptive', False)
+        accuracy          = kwargs.pop('accuracy', 1.0e-4)
+        maxiter           = kwargs.pop('maxiter', 20)
+        adaptive          = kwargs.pop('adaptive', False)
+        detect_divergence = kwargs.pop('detect_divergence', True)
+        quiet             = kwargs.pop('quiet', False)
 
         # initialize iterative process:
-        x0, y0 = self.wcs_sky2pix(ra, dec, origin) # initial approximation (WCS based only)
+        x0, y0 = self.wcs_sky2pix(ra, dec, origin) # <-- initial approximation
+                                                  #     (WCS based only)
 
         # see if iterative solution is required (when any of the
         # non-CD-matrix corrections are present). If not required
@@ -585,7 +720,8 @@ accuracy after 3 iterations.
         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:
-            # no non-WCS corrections are detected - return initial approximation
+            # no non-WCS corrections are detected - return
+            # initial approximation
             if vect1D:
                 return [x0, y0]
             else:
@@ -608,8 +744,8 @@ accuracy after 3 iterations.
         y -= dy
 
         # norn (L2) squared of the correction:
-        dn2prev = dx**2+dy**2
-        dn2     = dn2prev
+        dn2prev   = dx**2+dy**2
+        dn2       = dn2prev
 
         # process all coordinates simultaneously:
         iterlist  = range(1, maxiter+1)
@@ -617,7 +753,12 @@ accuracy after 3 iterations.
         ind       = None
         inddiv    = None
 
-        divergent = False
+        npts      = x.shape[0]
+        
+        # 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')
 
         if not adaptive:
             for k in iterlist:
@@ -625,66 +766,113 @@ accuracy after 3 iterations.
                 if np.max(dn2) < accuracy2:
                     break
 
-                # check for divergence:
-                inddiv, = np.where((dn2 > dn2prev) & (dn2 >= accuracy2))
-                if inddiv.shape[0] > 0:
-                    divergent = True
-                    break
-
                 # find correction to the previous solution:
                 dx, dy = self.pix2foc(x, y, origin)
                 # If pix2foc does not apply all the required distortion
                 # corrections then replace the above line with:
                 #r0, d0 = self.all_pix2sky(x, y, origin)
                 #dx, dy = self.wcs_sky2pix(r0, d0, origin )
-
                 dx -= x0
                 dy -= y0
 
+                # update norn (L2) squared of the correction:
+                dn2 = dx**2+dy**2
+
+                # check for divergence (we do this in two stages
+                # to optimize performance for the most common
+                # scenario when succesive approximations converge):
+                if detect_divergence:
+                    ind, = np.where(dn2 <= dn2prev)
+                    if ind.shape[0] < npts:
+                        inddiv, = np.where(
+                            np.logical_and(dn2 > dn2prev, dn2 >= accuracy2))
+                        if inddiv.shape[0] > 0:
+                            # apply correction only to the converging points:
+                            x[ind] -= dx[ind]
+                            y[ind] -= dy[ind]
+                            # switch to adaptive iterations:
+                            ind, = np.where((dn2 >= accuracy2) & \
+                                        (dn2 <= dn2prev) & np.isfinite(dn2))
+                            iterlist = iterlist[k:]
+                            adaptive = True
+                            break
+                    #dn2prev[ind] = dn2[ind]
+                    dn2prev = dn2
+
                 # apply correction:
                 x -= dx
                 y -= dy
 
-                # update norn (L2) squared of the correction:
-                dn2prev = dn2.copy()
-                dn2     = dx**2+dy**2
+        if adaptive:
 
-        else:
-            ind, = np.where(dn2 >= accuracy2)
+            if ind is None:
+                ind = np.asarray(range(npts), dtype=np.int64)
 
             for k in iterlist:
                 # check convergence:
                 if ind.shape[0] == 0:
                     break
 
-                # check for divergence:
-                inddiv = ind[np.where(dn2[ind] > dn2prev[ind])]
-                if inddiv.shape[0] > 0:
-                    divergent = True
-                    break
-
                 # find correction to the previous solution:
                 dx[ind], dy[ind] = self.pix2foc(x[ind], y[ind], origin)
                 # If pix2foc does not apply all the required distortion
                 # corrections then replace the above line with:
                 #r0[ind], d0[ind] = self.all_pix2sky(x[ind], y[ind], origin)
-                #dx[ind], dy[ind] = self.wcs_sky2pix(r0[ind], d0[ind], origin )
+                #dx[ind], dy[ind] = self.wcs_sky2pix(r0[ind], d0[ind], origin)
                 dx[ind] -= x0[ind]
                 dy[ind] -= y0[ind]
 
+                # update norn (L2) squared of the correction:
+                dn2 = dx**2+dy**2
+
+                # update indices of elements that still need correction:
+                if detect_divergence:
+                    ind, = np.where((dn2 >= accuracy2) & (dn2 <= dn2prev))
+                    dn2prev[ind] = dn2[ind]
+                else:
+                    ind, = np.where(dn2 >= accuracy2)
+                    #ind = ind[np.where(dx[ind]**2+dy[ind]**2 >= accuracy2)]
+
                 # apply correction:
                 x[ind] -= dx[ind]
                 y[ind] -= dy[ind]
 
-                # update norn (L2) squared of the correction:
-                dn2prev  = dn2.copy()
-                dn2 = dx**2+dy**2
+        #FINAL DETECTION OF INVALID, DIVERGING, AND FAILED-TO-CONVERGE 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:
+                ind = None
 
-                # update indices of elements that still need correction:
-                ind, = np.where(dn2 >= accuracy2)
-                #ind = ind[np.where(dx[ind]**2+dy[ind]**2 >= accuracy2)]
+        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
 
-        if k >= maxiter and not quiet:
+        if (ind is not None or inddiv is not None) and not quiet:
             if vect1D:
                 sol  = [x, y]
                 err  = [np.abs(dx), np.abs(dy)]
@@ -692,31 +880,31 @@ accuracy after 3 iterations.
                 sol  = np.dstack( [x, y] )[0]
                 err  = np.dstack( [np.abs(dx), np.abs(dy)] )[0]
 
-            if ind is None:
-                ind, = np.where(dx**2+dy**2 >= accuracy2)
-
-            if inddiv is None:
-                inddiv, = np.where(dn2[ind] > dn2prev[ind])
-
-            if ind.shape[0] == 0:
-                ind    = None
-                inddiv = None
-
-            elif inddiv.shape[0] == 0:
-                inddiv = None
-
-            assert(ind is not None or inddiv is not None) # <-- sanity check
+            # restore previous numpy error settings:
+            np.seterr(invalid = old_invalid, over = old_over)
+            
+            if inddiv is None:                
+                raise NoConvergence("'HSTWCS.all_sky2pix' failed to "        \
+                    "converge to the requested accuracy after {:d} "         \
+                    "iterations.".format(k), best_solution = sol,            \
+                    accuracy = err, niter = k, failed2converge = ind,        \
+                    divergent = None)
+            else:
+                raise NoConvergence("'HSTWCS.all_sky2pix' failed to "        \
+                    "converge to the requested accuracy.{0:s}"               \
+                    "After {1:d} iterations, the solution is diverging "     \
+                    "at least for one input point."                          \
+                    .format(os.linesep, k), best_solution = sol,             \
+                    accuracy = err, niter = k, failed2converge = ind,        \
+                    divergent = inddiv)
 
-            raise NoConvergence("'HSTWCS.all_sky2pix' failed to converge to "\
-                                "requested accuracy after {:d} iterations."  \
-                                .format(k), best_solution = sol,             \
-                                accuracy = err, niter = k,                   \
-                                nonconvergent = ind, divergent = inddiv)
+        # restore previous numpy error settings:
+        np.seterr(invalid = old_invalid, over = old_over)
 
         if vect1D:
             return [x, y]
         else:
-            return np.dstack([x,y])[0]
+            return np.dstack( [x, y] )[0]
 
 
     def _updatehdr(self, ext_hdr):