From 95776c98790b249e62f78216cdd8eff5b5c8cfca Mon Sep 17 00:00:00 2001 From: Debora Pelliccia Date: Wed, 15 May 2024 17:21:03 -1000 Subject: [PATCH 1/2] add kind=stable to argsort (cherry picked from commit e6abfcd5e66d9777c680fd22165721c963ba6dbd) --- pypeit/bitmask.py | 2 +- pypeit/bspline/bspline.py | 6 +++--- pypeit/coadd3d.py | 4 ++-- pypeit/core/arc.py | 4 ++-- pypeit/core/datacube.py | 4 ++-- pypeit/core/extract.py | 18 +++++++++--------- pypeit/core/findobj_skymask.py | 14 +++++++------- pypeit/core/fitting.py | 2 +- pypeit/core/flat.py | 6 +++--- pypeit/core/pydl.py | 6 +++--- pypeit/core/skysub.py | 12 ++++++------ pypeit/core/slitdesign_matching.py | 4 ++-- pypeit/core/telluric.py | 2 +- pypeit/core/wavecal/autoid.py | 2 +- pypeit/edgetrace.py | 4 ++-- pypeit/flatfield.py | 8 ++++---- pypeit/sampling.py | 2 +- pypeit/utils.py | 4 ++-- 18 files changed, 52 insertions(+), 52 deletions(-) diff --git a/pypeit/bitmask.py b/pypeit/bitmask.py index fd35e61c0f..cfa4eae9ca 100644 --- a/pypeit/bitmask.py +++ b/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]) diff --git a/pypeit/bspline/bspline.py b/pypeit/bspline/bspline.py index d4efc6c10b..2867460710 100644 --- a/pypeit/bspline/bspline.py +++ b/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. diff --git a/pypeit/coadd3d.py b/pypeit/coadd3d.py index 225747d76d..5638cb6f77 100644 --- a/pypeit/coadd3d.py +++ b/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) diff --git a/pypeit/core/arc.py b/pypeit/core/arc.py index e69230cf28..086f7deb04 100644 --- a/pypeit/core/arc.py +++ b/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 diff --git a/pypeit/core/datacube.py b/pypeit/core/datacube.py index 3f9625ba9b..9288e06695 100644 --- a/pypeit/core/datacube.py +++ b/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 diff --git a/pypeit/core/extract.py b/pypeit/core/extract.py index 1c3777d61b..a820b27be5 100644 --- a/pypeit/core/extract.py +++ b/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]) diff --git a/pypeit/core/findobj_skymask.py b/pypeit/core/findobj_skymask.py index fe56f70bc4..24979b84cc 100644 --- a/pypeit/core/findobj_skymask.py +++ b/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 diff --git a/pypeit/core/fitting.py b/pypeit/core/fitting.py index 7f4346ebe4..eddc49ec07 100644 --- a/pypeit/core/fitting.py +++ b/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'] diff --git a/pypeit/core/flat.py b/pypeit/core/flat.py index 46495f56f2..bac66b4c9a 100644 --- a/pypeit/core/flat.py +++ b/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) diff --git a/pypeit/core/pydl.py b/pypeit/core/pydl.py index b7916dfbd4..97c093ccea 100644 --- a/pypeit/core/pydl.py +++ b/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 # diff --git a/pypeit/core/skysub.py b/pypeit/core/skysub.py index a5114d886c..d92a734b82 100644 --- a/pypeit/core/skysub.py +++ b/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) diff --git a/pypeit/core/slitdesign_matching.py b/pypeit/core/slitdesign_matching.py index eda3db0682..f12ce01b45 100644 --- a/pypeit/core/slitdesign_matching.py +++ b/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 Date: Tue, 13 Aug 2024 13:39:20 -1000 Subject: [PATCH 2/2] added to release --- doc/releases/1.16.1dev.rst | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/doc/releases/1.16.1dev.rst b/doc/releases/1.16.1dev.rst index a223665b4a..a725a302c7 100644 --- a/doc/releases/1.16.1dev.rst +++ b/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. \ No newline at end of file +- 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. + + +