Merge pull request #1840 from pypeit/argsort_fix

Numpy Argsort fix

doc/releases/1.16.1dev.rst

@@ -55,4 +55,10 @@ Bug Fixes
 - Fixed a fault caused when all frames in a pypeit file are identified as being
   part of ``all`` calibration groups.
 - Allow for empty 2D wavecal solution in HDU extension of WaveCalib file
-- Fixed a bug in the ginga display function, when the user doesn't provide the `trc_name` argument.
+- Fixed a bug in the ginga display function, when the user doesn't provide the `trc_name` argument.
+- Fix a MAJOR BUT SUBTLE bug in the use of ``numpy.argsort``. When using ``numpy.argsort``
+  the parameter kind='stable' should be used to ensure that a sorting algorithm more robust
+  than "quicksort" is used.
+
+
+

pypeit/bitmask.py

@@ -575,7 +575,7 @@ def from_header(cls, hdr, prefix=None):
         # Fill in any missing bits
         keys, values, descr = cls._fill_sequence(keys, values, descr=descr)
         # Make sure the bits are sorted
-        srt = numpy.argsort(values)
+        srt = numpy.argsort(values, kind='stable')
         # Instantiate the BitMask
         return cls(keys[srt], descr=descr[srt])
 

pypeit/bspline/bspline.py

@@ -532,7 +532,7 @@ def value(self, x, x2=None, action=None, lower=None, upper=None):
             is good).
         """
         # TODO: Is the sorting necessary?
-        xsort = x.argsort()
+        xsort = x.argsort(kind='stable')
         if action is None:
             action, lower, upper = self.action(x[xsort], x2=None if x2 is None else x2[xsort])
         else:
@@ -550,12 +550,12 @@ def value(self, x, x2=None, action=None, lower=None, upper=None):
         mask[(x < gb[self.nord-1]) | (x > gb[n])] = False
         hmm = (np.diff(goodbk) > 2).nonzero()[0]
         if hmm.size == 0:
-            return yfit[np.argsort(xsort)], mask
+            return yfit[np.argsort(xsort, kind='stable')], mask
 
         for jj in range(hmm.size):
             mask[(x >= self.breakpoints[goodbk[hmm[jj]]])
                     & (x <= self.breakpoints[goodbk[hmm[jj]+1]-1])] = False
-        return yfit[np.argsort(xsort)], mask
+        return yfit[np.argsort(xsort, kind='stable')], mask
 
     def maskpoints(self, err):
         """Perform simple logic of which breakpoints to mask.

pypeit/coadd3d.py

@@ -1116,11 +1116,11 @@ def load(self):
             dwav_ext = dwaveimg[onslit_gpm]
 
             # For now, work in sorted wavelengths
-            wvsrt = np.argsort(wave_ext)
+            wvsrt = np.argsort(wave_ext, kind='stable')
             wave_sort = wave_ext[wvsrt]
             dwav_sort = dwav_ext[wvsrt]
             # Here's an array to get back to the original ordering
-            resrt = np.argsort(wvsrt)
+            resrt = np.argsort(wvsrt, kind='stable')
 
             # Compute the DAR correction
             cosdec = np.cos(self.ifu_dec[ff] * np.pi / 180.0)

pypeit/core/arc.py

@@ -644,7 +644,7 @@ def detect_peaks(x, mph=None, mpd=1, threshold=0, edge='rising',
         ind = np.delete(ind, np.where(dx < threshold)[0])
     # detect small peaks closer than minimum peak distance
     if ind.size and mpd > 1:
-        ind = ind[np.argsort(x[ind])][::-1]  # sort ind by peak height
+        ind = ind[np.argsort(x[ind], kind='stable')][::-1]  # sort ind by peak height
         idel = np.zeros(ind.size, dtype=bool)
         for i in range(ind.size):
             if not idel[i]:
@@ -808,7 +808,7 @@ def iter_continuum(spec, gpm=None, fwhm=4.0, sigthresh = 2.0, sigrej=3.0, niter_
             #cont_mask = np.ones_like(cont_mask) & gpm
             # Short circuit the masking and just mask the 0.70 most offending pixels
             peak_mask_ind = np.where(np.logical_not(peak_mask) & gpm)[0]
-            isort = np.argsort(np.abs(spec[peak_mask_ind]))[::-1]
+            isort = np.argsort(np.abs(spec[peak_mask_ind]), kind='stable')[::-1]
             peak_mask_new = np.ones_like(peak_mask)
             peak_mask_new[peak_mask_ind[isort[0:max_nmask]]] = False
             cont_mask = peak_mask_new & gpm

pypeit/core/datacube.py

@@ -1754,7 +1754,7 @@ def subpixellate(output_wcs, bins, sciImg, ivarImg, waveImg, slitid_img_gpm, wgh
             yspl = this_tilts[wpix] * (this_slits.nspec - 1)
             tiltpos = np.add.outer(yspl, spec_y).flatten()
             wspl = this_wav[this_sl]
-            asrt = np.argsort(yspl)
+            asrt = np.argsort(yspl, kind='stable')
             wave_spl = interp1d(yspl[asrt], wspl[asrt], kind='linear', bounds_error=False, fill_value='extrapolate')
             # Calculate the wavelength at each subpixel
             this_wave_subpix = wave_spl(tiltpos)
@@ -1769,7 +1769,7 @@ def subpixellate(output_wcs, bins, sciImg, ivarImg, waveImg, slitid_img_gpm, wgh
             # Transform this to spatial location
             spatpos_subpix = _astrom_trans[fr].transform(sl, spat_xx, spec_yy)
             spatpos = _astrom_trans[fr].transform(sl, wpix[1], wpix[0])
-            ssrt = np.argsort(spatpos)
+            ssrt = np.argsort(spatpos, kind='stable')
             # Initialize the voxel coordinates for each spec2D pixel
             vox_coord = np.full((numpix, num_all_subpixels, 3), -1, dtype=float)
             # Loop over the subslices

pypeit/core/extract.py

@@ -644,7 +644,7 @@ def qa_fit_profile(x_tot, y_tot, model_tot, l_limit = None,
             closest = (dist[close]).argmin()
             model_samp[i] = (model[close])[closest]
         if nclose > 3:
-            s = yclose.argsort()
+            s = yclose.argsort(kind='stable')
             y50[i] = yclose[s[int(np.rint((nclose - 1)*0.5))]]
             y80[i] = yclose[s[int(np.rint((nclose - 1)*0.8))]]
             y20[i] = yclose[s[int(np.rint((nclose - 1)*0.2))]]
@@ -656,7 +656,7 @@ def qa_fit_profile(x_tot, y_tot, model_tot, l_limit = None,
     else:
         ax.plot(plot_mid, y50, marker='o', color='lime', markersize=2, fillstyle = 'full', linestyle='None')
 
-    isort = x.argsort()
+    isort = x.argsort(kind='stable')
     ax.plot(x[isort], model[isort], color='red', linewidth=1.0)
 
 
@@ -905,7 +905,7 @@ def fit_profile(image, ivar, waveimg, thismask, spat_img, trace_in, wave,
     spline_flux1[ispline] = spline_tmp
     cont_tmp, _ = c_answer.value(wave[ispline])
     cont_flux1[ispline] = cont_tmp
-    isrt = np.argsort(wave[indsp])
+    isrt = np.argsort(wave[indsp], kind='stable')
     s2_1_interp = scipy.interpolate.interp1d(wave[indsp][isrt], sn2[isrt],assume_sorted=False, bounds_error=False,fill_value = 0.0)
     sn2_1[ispline] = s2_1_interp(wave[ispline])
     bmask = np.zeros(nspec,dtype='bool')
@@ -952,7 +952,7 @@ def fit_profile(image, ivar, waveimg, thismask, spat_img, trace_in, wave,
         # Create the normalized object image
         if np.any(totmask):
             igd = (wave >= wave_min) & (wave <= wave_max)
-            isrt1 = np.argsort(wave[igd])
+            isrt1 = np.argsort(wave[igd], kind='stable')
             #plt.plot(wave[igd][isrt1], spline_flux1[igd][isrt1])
             #plt.show()
             spline_img_interp = scipy.interpolate.interp1d(wave[igd][isrt1],spline_flux1[igd][isrt1],assume_sorted=False,
@@ -1018,7 +1018,7 @@ def fit_profile(image, ivar, waveimg, thismask, spat_img, trace_in, wave,
         inside, = np.where(IN_PIX.flatten())
 
 
-    si = inside[np.argsort(sigma_x.flat[inside])]
+    si = inside[np.argsort(sigma_x.flat[inside], kind='stable')]
     sr = si[::-1]
 
     bset, bmask = fitting.iterfit(sigma_x.flat[si],norm_obj.flat[si], invvar = norm_ivar.flat[si],
@@ -1075,7 +1075,7 @@ def fit_profile(image, ivar, waveimg, thismask, spat_img, trace_in, wave,
         return (profile_model, trace_in, fwhmfit, med_sn2)
 
     sigma_iter = 3
-    isort = (xtemp.flat[si[inside]]).argsort()
+    isort = (xtemp.flat[si[inside]]).argsort(kind='stable')
     inside = si[inside[isort]]
     pb = np.ones(inside.size)
 
@@ -1152,7 +1152,7 @@ def fit_profile(image, ivar, waveimg, thismask, spat_img, trace_in, wave,
 
         # Update the profile B-spline fit for the next iteration
         if iiter < sigma_iter-1:
-            ss = sigma_x.flat[inside].argsort()
+            ss = sigma_x.flat[inside].argsort(kind='stable')
             pb = (np.outer(area, np.ones(nspat,dtype=float))).flat[inside]
             keep = (bkpt >= sigma_x.flat[inside].min()) & (bkpt <= sigma_x.flat[inside].max())
             if keep.sum() == 0:
@@ -1176,7 +1176,7 @@ def fit_profile(image, ivar, waveimg, thismask, spat_img, trace_in, wave,
         xnew = trace_in
 
     fwhmfit = sigma*2.3548
-    ss=sigma_x.flatten().argsort()
+    ss=sigma_x.flatten().argsort(kind='stable')
     inside, = np.where((sigma_x.flat[ss] >= min_sigma) &
                        (sigma_x.flat[ss] <= max_sigma) &
                        mask[ss] &
@@ -1195,7 +1195,7 @@ def fit_profile(image, ivar, waveimg, thismask, spat_img, trace_in, wave,
     sigma_x_igood = sigma_x.flat[igood]
     yfit_out, _  = bset.value(sigma_x_igood)
     full_bsp[igood] = yfit_out
-    isrt2 = sigma_x_igood.argsort()
+    isrt2 = sigma_x_igood.argsort(kind='stable')
     (peak, peak_x, lwhm, rwhm) = findfwhm(yfit_out[isrt2] - median_fit, sigma_x_igood[isrt2])
 
 

pypeit/core/findobj_skymask.py

@@ -523,7 +523,7 @@ def ech_fill_in_orders(sobjs:specobjs.SpecObjs,
     uni_frac = gfrac[uni_ind]
 
     # Sort with respect to fractional slit location to guarantee that we have a similarly sorted list of objects later
-    isort_frac = uni_frac.argsort()
+    isort_frac = uni_frac.argsort(kind='stable')
     uni_obj_id = uni_obj_id[isort_frac]
     uni_frac = uni_frac[isort_frac]
 
@@ -788,7 +788,7 @@ def ech_cutobj_on_snr(
 
     ## Loop over objects from highest SNR to lowest SNR. Apply the S/N constraints. Once we hit the maximum number
     # objects requested exit, except keep the hand apertures that were requested.
-    isort_SNR_max = np.argsort(np.median(SNR_arr,axis=0))[::-1]
+    isort_SNR_max = np.argsort(np.median(SNR_arr,axis=0), kind='stable')[::-1]
     for iobj in isort_SNR_max:
         hand_ap_flag = hand_flag[iobj]
         SNR_constraint = (SNR_arr[:,iobj].max() > max_snr) or (
@@ -800,7 +800,7 @@ def ech_cutobj_on_snr(
             sobjs_keep = sobjs_align[ikeep].copy()
             sobjs_keep.ECH_OBJID = iobj_keep
             sobjs_keep.OBJID = iobj_keep
-            sobjs_trim.add_sobj(sobjs_keep[np.argsort(sobjs_keep.SLITID)])
+            sobjs_trim.add_sobj(sobjs_keep[np.argsort(sobjs_keep.SLITID, kind='stable')])
             iobj_keep += 1
             if not hand_ap_flag:
                 iobj_keep_not_hand += 1
@@ -1700,7 +1700,7 @@ def orig_ech_objfind(image, ivar, slitmask, slit_left, slit_righ, order_vec, mas
     uni_frac = gfrac[uni_ind]
 
     # Sort with respect to fractional slit location to guarantee that we have a similarly sorted list of objects later
-    isort_frac = uni_frac.argsort()
+    isort_frac = uni_frac.argsort(kind='stable')
     uni_obj_id = uni_obj_id[isort_frac]
     uni_frac = uni_frac[isort_frac]
 
@@ -1865,7 +1865,7 @@ def orig_ech_objfind(image, ivar, slitmask, slit_left, slit_righ, order_vec, mas
 
     ## Loop over objects from highest SNR to lowest SNR. Apply the S/N constraints. Once we hit the maximum number
     # objects requested exit, except keep the hand apertures that were requested.
-    isort_SNR_max = np.argsort(np.median(SNR_arr,axis=0))[::-1]
+    isort_SNR_max = np.argsort(np.median(SNR_arr,axis=0), kind='stable')[::-1]
     for iobj in isort_SNR_max:
         hand_ap_flag = int(np.round(slitfracpos_arr[0, iobj]*1000)) in hand_frac
         SNR_constraint = (SNR_arr[:,iobj].max() > max_snr) or (np.sum(SNR_arr[:,iobj] > min_snr) >= nabove_min_snr)
@@ -1879,7 +1879,7 @@ def orig_ech_objfind(image, ivar, slitmask, slit_left, slit_righ, order_vec, mas
 #            for spec in sobjs_keep:
 #                spec.ECH_OBJID = iobj_keep
 #                #spec.OBJID = iobj_keep
-            sobjs_trim.add_sobj(sobjs_keep[np.argsort(sobjs_keep.ECH_ORDERINDX)])
+            sobjs_trim.add_sobj(sobjs_keep[np.argsort(sobjs_keep.ECH_ORDERINDX, kind='stable')])
             iobj_keep += 1
             if not hand_ap_flag:
                 iobj_keep_not_hand += 1
@@ -2684,7 +2684,7 @@ def objs_in_slit(image, ivar, thismask, slit_left, slit_righ,
     # Sort objects according to their spatial location
     nobj = len(sobjs)
     spat_pixpos = sobjs.SPAT_PIXPOS
-    sobjs = sobjs[spat_pixpos.argsort()]
+    sobjs = sobjs[spat_pixpos.argsort(kind='stable')]
     # Assign integer objids
     sobjs.OBJID = np.arange(nobj) + 1
 

pypeit/core/fitting.py

@@ -943,7 +943,7 @@ def iterfit(xdata, ydata, invvar=None, inmask=None, upper=5, lower=5, x2=None,
         outmask = True
     else:
         outmask = np.ones(invvar.shape, dtype='bool')
-    xsort = xdata.argsort()
+    xsort = xdata.argsort(kind='stable')
     maskwork = (outmask & inmask & (invvar > 0.0))[xsort]  # `maskwork` is in xsort order
     if 'oldset' in kwargs_bspline:
         sset = kwargs_bspline['oldset']

pypeit/core/flat.py

@@ -67,7 +67,7 @@ def sorted_flat_data(data, coo, gpm=None):
     # np.argsort sorts the data over the last axis. To avoid coo[gpm]
     # returning an array (which will happen if the gpm is not provided
     # as an argument), all the arrays are explicitly flattened.
-    srt = np.argsort(coo[gpm].ravel())
+    srt = np.argsort(coo[gpm].ravel(), kind='stable')
     coo_data = coo[gpm].ravel()[srt]
     flat_data = data[gpm].ravel()[srt]
     return gpm, srt, coo_data, flat_data
@@ -241,7 +241,7 @@ def construct_illum_profile(norm_spec, spat_coo, slitwidth, spat_gpm=None, spat_
         plt.show()
 
     # Include the rejected data in the full image good-pixel mask
-    _spat_gpm[_spat_gpm] = spat_gpm_data_raw[np.argsort(spat_srt)]
+    _spat_gpm[_spat_gpm] = spat_gpm_data_raw[np.argsort(spat_srt, kind='stable')]
     # Recreate the illumination profile data
     _spat_gpm, spat_srt, spat_coo_data, spat_flat_data_raw \
             = sorted_flat_data(norm_spec, spat_coo, gpm=_spat_gpm)
@@ -434,7 +434,7 @@ def poly_map(rawimg, rawivar, waveimg, slitmask, slitmask_trim, modelimg, deg=3,
     slitmask_spatid = np.sort(slitmask_spatid[slitmask_spatid > 0])
 
     # Create a spline between the raw data and the error
-    flxsrt = np.argsort(np.ravel(rawimg))
+    flxsrt = np.argsort(np.ravel(rawimg), kind='stable')
     spl = scipy.interpolate.interp1d(np.ravel(rawimg)[flxsrt], np.ravel(rawivar)[flxsrt], kind='linear',
                                      bounds_error=False, fill_value=0.0, assume_sorted=True)
     modelmap = np.ones_like(rawimg)

pypeit/core/pydl.py

@@ -55,7 +55,7 @@ def djs_maskinterp1(yval, mask, xval=None, const=False):
             if igood[ngood-1] != ny-1:
                 ynew[igood[ngood-1]+1:ny] = ynew[igood[ngood-1]]
     else:
-        ii = xval.argsort()
+        ii = xval.argsort(kind='stable')
         ibad = (mask[ii] != 0).nonzero()[0]
         igood = (mask[ii] == 0).nonzero()[0]
         ynew[ii[ibad]] = np.interp(xval[ii[ibad]], xval[ii[igood]],
@@ -882,7 +882,7 @@ def djs_reject(data, model, outmask=None, inmask=None,
                         # Test if too many points rejected in this group.
                         #
                         if np.sum(badness[jj] != 0) > maxrej1[iloop]:
-                            isort = badness[jj].argsort()
+                            isort = badness[jj].argsort(kind='stable')
                             #
                             # Make the following points good again.
                             #
@@ -1656,7 +1656,7 @@ def spherematch(ra1, dec1, ra2, dec2, matchlength, chunksize=None,
     omatch1 = np.array(match1)
     omatch2 = np.array(match2)
     odistance12 = np.array(distance12)
-    s = odistance12.argsort()
+    s = odistance12.argsort(kind='stable')
     #
     # Retain only desired matches
     #

pypeit/core/skysub.py

@@ -147,7 +147,7 @@ def global_skysub(image, ivar, tilts, thismask, slit_left, slit_righ, inmask=Non
         return np.zeros(np.sum(thismask))
 
     # Sub arrays
-    isrt = np.argsort(piximg[thismask])
+    isrt = np.argsort(piximg[thismask], kind='stable')
     pix = piximg[thismask][isrt]
     sky = image[thismask][isrt]
     sky_ivar = ivar[thismask][isrt]
@@ -302,7 +302,7 @@ def skyoptimal(piximg, data, ivar, oprof, sigrej=3.0, npoly=1, spatial_img=None,
         whether a pixel is good (True) or was masked (False).
     """
 
-    sortpix = piximg.argsort()
+    sortpix = piximg.argsort(kind='stable')
 
     nx = data.size
     nc = oprof.shape[0]
@@ -328,7 +328,7 @@ def skyoptimal(piximg, data, ivar, oprof, sigrej=3.0, npoly=1, spatial_img=None,
     indx, = np.where(ivar[sortpix] > 0.0)
     ngood = indx.size
     good = sortpix[indx]
-    good = good[piximg[good].argsort()]
+    good = good[piximg[good].argsort(kind='stable')]
     relative, = np.where(relative_mask[good])
 
     gpm = np.zeros(piximg.shape, dtype=bool)
@@ -428,7 +428,7 @@ def optimal_bkpts(bkpts_optimal, bsp_min, piximg, sampmask, samp_frac=0.80,
     """
 
     pix = piximg[sampmask]
-    isrt = pix.argsort()
+    isrt = pix.argsort(kind='stable')
     pix = pix[isrt]
     piximg_min = pix.min()
     piximg_max = pix.max()
@@ -1319,11 +1319,11 @@ def ech_local_skysub_extract(sciimg, sciivar, fullmask, tilts, waveimg,
 
     # Compute the average SNR and find the brightest object
     snr_bar = np.mean(order_snr,axis=0)
-    srt_obj = snr_bar.argsort()[::-1]
+    srt_obj = snr_bar.argsort(kind='stable')[::-1]
     ibright = srt_obj[0] # index of the brightest object
 
     # Now extract the orders in descending order of S/N for the brightest object
-    srt_order_snr = order_snr[:,ibright].argsort()[::-1]
+    srt_order_snr = order_snr[:,ibright].argsort(kind='stable')[::-1]
     fwhm_here = np.zeros(norders)
     fwhm_was_fit = np.zeros(norders,dtype=bool)
 

pypeit/core/slitdesign_matching.py

@@ -74,7 +74,7 @@ def best_offset(x_det, x_model, step=1, xlag_range=None):
     for j in range(xlag.size):
         x_det_lag = x_det+xlag[j]
         join = np.ma.concatenate([x_det_lag, x_model_trim])
-        sind = np.argsort(join)
+        sind = np.argsort(join, kind='stable')
         nj = sind.size
         w1 = np.where(sind < x_det.size)
 
@@ -90,7 +90,7 @@ def best_offset(x_det, x_model, step=1, xlag_range=None):
             offs = np.amin(np.absolute(x_det_lag[:, None] - x_model_trim[None, :]), axis=1)
 
             # use only nbest best matches
-            soffs = np.argsort(np.abs(offs))
+            soffs = np.argsort(np.abs(offs), kind='stable')
             nbest2 = nbest if nbest<x_det.size else x_det.size
             offs = offs[soffs[0:nbest2]]
 

pypeit/core/telluric.py

@@ -2803,6 +2803,6 @@ def sort_telluric(self):
                 tell_med[iord] = np.mean(np.exp(-np.sinh(tell_model_mean)))
 
         # Perform fits in order of telluric strength
-        return tell_med.argsort()
+        return tell_med.argsort(kind='stable')
 
 

pypeit/core/wavecal/autoid.py

@@ -2339,7 +2339,7 @@ def cross_match(self, good_fit, detections):
                 wvc_gd_jfh[cntr] = wave_soln[self._npix//2]
                 dsp_gd_jfh[cntr]= np.median(wave_soln - np.roll(wave_soln,1))
                 cntr += 1
-        srt = np.argsort(wvc_gd_jfh)
+        srt = np.argsort(wvc_gd_jfh, kind='stable')
         sort_idx = idx_gd[srt]
         sort_wvc = wvc_gd[srt]
         sort_dsp = dsp_gd[srt]