From 0e87a0f12b1125af73ed91b02bfbfa5229159162 Mon Sep 17 00:00:00 2001 From: Niv Drory Date: Tue, 10 Sep 2024 10:04:59 -0500 Subject: [PATCH 1/5] replace numpy median/mean/... with faster bottleneck calls --- python/lvmdrp/core/cube.py | 2 +- python/lvmdrp/core/fiberrows.py | 3 +- python/lvmdrp/core/fluxcal.py | 8 +-- python/lvmdrp/core/image.py | 22 +++---- python/lvmdrp/core/rss.py | 8 +-- python/lvmdrp/core/spectrum1d.py | 30 +++++----- python/lvmdrp/functions/cubeMethod.py | 4 +- python/lvmdrp/functions/imageMethod.py | 74 ++++++++++++------------ python/lvmdrp/functions/plotMethod.py | 4 +- python/lvmdrp/functions/rssMethod.py | 46 +++++++-------- python/lvmdrp/functions/run_calseq.py | 2 +- python/lvmdrp/functions/run_twilights.py | 9 +-- 12 files changed, 107 insertions(+), 105 deletions(-) diff --git a/python/lvmdrp/core/cube.py b/python/lvmdrp/core/cube.py index cac20a42..6260090c 100644 --- a/python/lvmdrp/core/cube.py +++ b/python/lvmdrp/core/cube.py @@ -412,7 +412,7 @@ def collapseCube( mask = None elif mode == "median": if numpy.sum(select_wave) > 0: - image = numpy.median(data[select_wave, :, :], 0) + image = bn.median(data[select_wave, :, :], 0) mask = image == 0 else: image = None diff --git a/python/lvmdrp/core/fiberrows.py b/python/lvmdrp/core/fiberrows.py index 1533c45b..a309c2ec 100644 --- a/python/lvmdrp/core/fiberrows.py +++ b/python/lvmdrp/core/fiberrows.py @@ -4,6 +4,7 @@ from tqdm import tqdm from copy import deepcopy as copy +import bottleneck as bn from lvmdrp import log from scipy import optimize from astropy.table import Table @@ -1268,7 +1269,7 @@ def smoothTraceDist( # offset1 = self._data[150, select_wave] - new_trace[150, select_wave] # offset2 = self._data[200, select_wave] - new_trace[200, select_wave] - offset_mean = numpy.median( + offset_mean = bn.median( self._data[:, select_wave] - new_trace[:, select_wave], axis=0 ) # computes that absolut trace position between the initially measured and estimated trace to compute the zero-point # offset_rms = numpy.std( diff --git a/python/lvmdrp/core/fluxcal.py b/python/lvmdrp/core/fluxcal.py index a5430d90..905288eb 100644 --- a/python/lvmdrp/core/fluxcal.py +++ b/python/lvmdrp/core/fluxcal.py @@ -12,7 +12,7 @@ from scipy import interpolate import requests import pandas as pd - +import bottleneck as bn import os.path as path import pathlib @@ -154,8 +154,8 @@ def mean_absolute_deviation(vals): Robust estimate of RMS - see https://en.wikipedia.org/wiki/Median_absolute_deviation """ - mval = np.nanmedian(vals) - rms = 1.4826 * np.nanmedian(np.abs(vals - mval)) + mval = bn.nanmedian(vals) + rms = 1.4826 * bn.nanmedian(np.abs(vals - mval)) return mval, rms # ok=np.abs(vals-mval)<4*rms @@ -312,7 +312,7 @@ def sky_flux_in_filter(cam, skyfibs, obswave, percentile=75): limidx = int(nfiber*percentile/100.0) skies = np.argsort(flux)[1:limidx] - return np.nanmedian(flux[skies]) + return bn.nanmedian(flux[skies]) def interpolate_mask(x, y, mask, kind="linear", fill_value=0): diff --git a/python/lvmdrp/core/image.py b/python/lvmdrp/core/image.py index 38ec322e..8c9dbd10 100644 --- a/python/lvmdrp/core/image.py +++ b/python/lvmdrp/core/image.py @@ -110,7 +110,7 @@ def _model_overscan(os_quad, axis=1, overscan_stat="biweight", threshold=None, m if overscan_stat == "biweight": stat = partial(biweight_location, ignore_nan=True) elif overscan_stat == "median": - stat = numpy.nanmedian + stat = bn.nanmedian else: warnings.warn( f"overscan statistic '{overscan_stat}' not implemented, " @@ -742,15 +742,15 @@ def measure_fiber_shifts(self, ref_image, columns=[500, 1000, 1500, 2000, 2500, shifts = numpy.zeros(len(columns)) for j,c in enumerate(columns): - s1 = numpy.nanmedian(ref_data[50:-50,c-column_width:c+column_width], axis=1) - s2 = numpy.nanmedian(self._data[50:-50,c-column_width:c+column_width], axis=1) - snr = numpy.sqrt(numpy.nanmedian(self._data[50:-50,c-column_width:c+column_width], axis=1)) + s1 = bn.nanmedian(ref_data[50:-50,c-column_width:c+column_width], axis=1) + s2 = bn.nanmedian(self._data[50:-50,c-column_width:c+column_width], axis=1) + snr = numpy.sqrt(bn.nanmedian(self._data[50:-50,c-column_width:c+column_width], axis=1)) min_snr = 5.0 - if numpy.nanmedian(snr) > min_snr: + if bn.nanmedian(snr) > min_snr: _, shifts[j], _ = _cross_match_float(s1, s2, numpy.array([1.0]), shift_range, gauss_window=[-3,3], min_peak_dist=5.0, ax=axs[j]) else: - comstr = f"low SNR (<={min_snr}) for thermal shift at column {c}: {numpy.nanmedian(snr):.4f}, assuming = 0.0" + comstr = f"low SNR (<={min_snr}) for thermal shift at column {c}: {bn.nanmedian(snr):.4f}, assuming = 0.0" log.warning(comstr) self.add_header_comment(comstr) shifts[j] = 0.0 @@ -2068,9 +2068,9 @@ def fitSpline(self, axis="y", degree=3, smoothing=0, use_weights=False, clip=Non # collapse groups into single pixel new_masked_pixels, new_data, new_vars = [], [], [] for group in groups: - new_masked_pixels.append(numpy.nanmean(masked_pixels[group])) - new_data.append(numpy.nanmedian(data[group])) - new_vars.append(numpy.nanmean(vars[group])) + new_masked_pixels.append(bn.nanmean(masked_pixels[group])) + new_data.append(bn.nanmedian(data[group])) + new_vars.append(bn.nanmean(vars[group])) masked_pixels = numpy.asarray(new_masked_pixels) data = numpy.asarray(new_data) vars = numpy.asarray(new_vars) @@ -2157,7 +2157,7 @@ def match_reference_column(self, ref_column=2000, ref_centroids=None, stretch_ra profile._data = numpy.nan_to_num(profile._data, nan=0, neginf=0, posinf=0) pixels = profile._pixels pixels = numpy.arange(pixels.size) - guess_heights = numpy.ones_like(ref_centroids) * numpy.nanmax(profile._data) + guess_heights = numpy.ones_like(ref_centroids) * bn.nanmax(profile._data) ref_profile = _spec_from_lines(ref_centroids, sigma=1.2, wavelength=pixels, heights=guess_heights) log.info(f"correcting guess positions for column {ref_column}") cc, bhat, mhat = _cross_match( @@ -3057,7 +3057,7 @@ def reject_cosmics(self, sigma_det=5, rlim=1.2, iterations=5, fwhm_gauss=[2.0,2. box_y = int(replace_box[1]) # define Laplacian convolution kernal - LA_kernel = numpy.array([[0, -1, 0], [-1, 4, -1], [0, -1, 0]])/4.0 + LA_kernel = 0.25*numpy.array([[0, -1, 0], [-1, 4, -1], [0, -1, 0]], dtype=numpy.float32) # Initiate image instances img_original = Image(data=self._data) diff --git a/python/lvmdrp/core/rss.py b/python/lvmdrp/core/rss.py index 3185827f..e34136a2 100644 --- a/python/lvmdrp/core/rss.py +++ b/python/lvmdrp/core/rss.py @@ -475,7 +475,7 @@ def from_channels(cls, rss_b, rss_r, rss_z, use_weights=True): new_data = bn.nansum(fluxes * weights, axis=0) new_lsf = bn.nansum(lsfs * weights, axis=0) new_error = numpy.sqrt(bn.nansum(vars, axis=0)) - new_mask = (numpy.nansum(masks, axis=0)>0) + new_mask = (bn.nansum(masks, axis=0)>0) if rss._sky is not None: new_sky = bn.nansum(skies * weights, axis=0) else: @@ -505,7 +505,7 @@ def from_channels(cls, rss_b, rss_r, rss_z, use_weights=True): new_data = bn.nanmean(fluxes, axis=0) new_lsf = bn.nanmean(lsfs, axis=0) new_error = numpy.sqrt(bn.nanmean(vars, axis=0)) - new_mask = numpy.nansum(masks, axis=0).astype("bool") + new_mask = bn.nansum(masks, axis=0).astype("bool") if skies.size != 0: new_sky = bn.nansum(skies, axis=0) else: @@ -1772,7 +1772,7 @@ def selectSpec(self, min=0, max=0, method="median"): if numpy.sum(goodpix) > 0: if method == "median": - collapsed[i] = numpy.median(spec._data[goodpix]) + collapsed[i] = bn.median(spec._data[goodpix]) elif method == "sum": collapsed[i] = numpy.sum(spec._data[goodpix]) elif method == "mean": @@ -3472,7 +3472,7 @@ def fit_field_gradient(self, wrange, poly_deg): x=fibermap["xpmm"].astype(float)[telescope=="Sci"] y=fibermap["ypmm"].astype(float)[telescope=="Sci"] - flux_med = numpy.nanmedian(flux) + flux_med = bn.nanmedian(flux) flux_fact = flux / flux_med select = numpy.isfinite(flux_fact) coeffs = polyfit2d(x[select], y[select], flux_fact[select], poly_deg) diff --git a/python/lvmdrp/core/spectrum1d.py b/python/lvmdrp/core/spectrum1d.py index 0c556a2a..4cebcb1b 100644 --- a/python/lvmdrp/core/spectrum1d.py +++ b/python/lvmdrp/core/spectrum1d.py @@ -38,7 +38,7 @@ def adaptive_smooth(data, start_width, end_width): start_index = max(0, i - half_width) end_index = min(n_points, i + half_width + 1) # Apply uniform filter to the local segment of the data - smoothed_data[i] = numpy.median(data[start_index:end_index]) + smoothed_data[i] = bn.median(data[start_index:end_index]) return smoothed_data def find_continuum(spec_s,niter=15,thresh=0.8,median_box_max=100,median_box_min=1): @@ -470,8 +470,8 @@ def wave_little_interpol(wavelist): # In overlap region patch in a linear scale with slightly different step. dw = overlap_end - overlap_start step = 0.5 * ( - numpy.mean(numpy.diff(wavelist[i])) - + numpy.mean(numpy.diff(wavelist[i + 1])) + bn.mean(numpy.diff(wavelist[i])) + + bn.mean(numpy.diff(wavelist[i + 1])) ) n_steps = int(dw / step + 0.5) @@ -2032,7 +2032,7 @@ def max(self): Pixel position of the maximum data value """ - max = numpy.nanmax(self._data) # get max + max = bn.nanmax(self._data) # get max select = self._data == max # select max value max_wave = self._wave[select][0] # get corresponding wavelength max_pos = self._pixels[select][0] # get corresponding position @@ -2054,7 +2054,7 @@ def min(self): Pixel position of the minimum data value """ - min = numpy.nanmin(self._data) # get min + min = bn.nanmin(self._data) # get min select = self._data == min # select min value min_wave = self._wave[select][0] # get corresponding waveength min_pos = self._pixels[select][0] # get corresponding position @@ -2104,7 +2104,7 @@ def resampleSpec( self._sky_error = numpy.flipud(self._sky_error) # case where input spectrum has more than half the pixels masked - if numpy.nansum(self._data) == 0.0 or ( + if bn.nansum(self._data) == 0.0 or ( self._mask is not None and numpy.sum(self._mask) > self._dim / 2 ): # all pixels masked @@ -3056,7 +3056,7 @@ def measureFWHMPeaks( pos_block = pos[ brackets[i] : brackets[i + 1] ] # cut out the corresponding peak positions - median_dist = numpy.nanmedian( + median_dist = bn.nanmedian( pos_block[1:] - pos_block[:-1] ) # compute median distance between peaks flux = ( @@ -3067,10 +3067,10 @@ def measureFWHMPeaks( ) # initial guess for the flux # compute lower and upper bounds of the positions for each block - lo = int(numpy.nanmin(pos_block) - median_dist) + lo = int(bn.nanmin(pos_block) - median_dist) if lo <= 0: lo = 0 - hi = int(numpy.nanmax(pos_block) + median_dist) + hi = int(bn.nanmax(pos_block) + median_dist) if hi >= self._wave[-1]: hi = self._wave[-1] @@ -3140,7 +3140,7 @@ def measureOffsetPeaks( pos_block = pos[blocks[i]] # cut out the corresponding peak positions pos_mask = good[blocks[i]] if numpy.sum(pos_mask) > 0: - median_dist = numpy.median( + median_dist = bn.median( pos_block[pos_mask][1:] - pos_block[pos_mask][:-1] ) # compute median distance between peaks flux = ( @@ -3188,7 +3188,7 @@ def measureOffsetPeaks( gaussians_offset.plot(self._wave[lo:hi], self._data[lo:hi]) offsets[i] = fit_par[-1] # get offset position - med_pos[i] = numpy.mean(self._wave[lo:hi]) + med_pos[i] = bn.mean(self._wave[lo:hi]) else: offsets[i] = 0.0 med_pos[i] = 0.0 @@ -3216,7 +3216,7 @@ def measureOffsetPeaks2( pos_mask = good[blocks[i]] pos_fwhm = fwhm[blocks[i]] if numpy.sum(pos_mask) > 0: - median_dist = numpy.median( + median_dist = bn.median( pos_block[pos_mask][1:] - pos_block[pos_mask][:-1] ) # compute median distance between peaks @@ -3252,7 +3252,7 @@ def measureOffsetPeaks2( max_flux[o] = numpy.sum(result[0]) find_max = numpy.argsort(max_flux)[-1] offsets[i] = offset[find_max] # get offset position - med_pos[i] = numpy.mean(self._wave[lo:hi]) + med_pos[i] = bn.mean(self._wave[lo:hi]) else: offsets[i] = 0.0 med_pos[i] = 0.0 @@ -3508,7 +3508,7 @@ def collapseSpec(self, method="mean", start=None, end=None, transmission_func=No if method != "mean" and method != "median" and method != "sum": raise ValueError("method must be either 'mean', 'median' or 'sum'") elif method == "mean": - flux = numpy.mean(self._data[select]) + flux = bn.mean(self._data[select]) if self._error is not None: error = numpy.sqrt( numpy.sum(self._error[select] ** 2) / numpy.sum(select) ** 2 @@ -3516,7 +3516,7 @@ def collapseSpec(self, method="mean", start=None, end=None, transmission_func=No else: error = None if self._sky is not None: - sky = numpy.mean(self._sky[select]) + sky = bn.mean(self._sky[select]) else: sky = None if self._sky_error is not None: diff --git a/python/lvmdrp/functions/cubeMethod.py b/python/lvmdrp/functions/cubeMethod.py index f0070e20..c5516276 100644 --- a/python/lvmdrp/functions/cubeMethod.py +++ b/python/lvmdrp/functions/cubeMethod.py @@ -1009,8 +1009,8 @@ def createTelluricCorrection_drp(cube_in, telluric_out, airmass, mask_telluric): line_par = stats.linregress( [mask_telluric[i * 2] - 10, mask_telluric[i * 2 + 1] + 10], [ - numpy.median(star_telluric1._data[select_blue]), - numpy.median(star_telluric1._data[select_red]), + bn.median(star_telluric1._data[select_blue]), + bn.median(star_telluric1._data[select_red]), ], ) star_telluric2._data[select_region] = ( diff --git a/python/lvmdrp/functions/imageMethod.py b/python/lvmdrp/functions/imageMethod.py index 8a39f186..5c190d92 100644 --- a/python/lvmdrp/functions/imageMethod.py +++ b/python/lvmdrp/functions/imageMethod.py @@ -131,8 +131,8 @@ def _nonlinearity_correction(ptc_params: None | numpy.ndarray, nominal_gain: flo gain_map = Image(data=1 / (a1 + a2*a3 * quadrant._data**(a3-1))) gain_map.setData(data=nominal_gain, select=numpy.isnan(gain_map._data), inplace=True) - gain_med = numpy.nanmedian(gain_map._data) - gain_min, gain_max = numpy.nanmin(gain_map._data), numpy.nanmax(gain_map._data) + gain_med = bn.nanmedian(gain_map._data) + gain_min, gain_max = bn.nanmin(gain_map._data), bn.nanmax(gain_map._data) log.info(f"gain map stats: {gain_med = :.2f} [{gain_min = :.2f}, {gain_max = :.2f}] ({nominal_gain = :.2f} e-/ADU)") else: log.warning("cannot apply non-linearity correction") @@ -391,8 +391,8 @@ def _fix_fiber_thermal_shifts(image, trace_cent, trace_width=None, trace_amp=Non # calculate thermal shifts column_shifts = image.measure_fiber_shifts(fiber_model, columns=columns, column_width=column_width, shift_range=shift_range, axs=axs_cc) # shifts stats - median_shift = numpy.nanmedian(column_shifts, axis=0) - std_shift = numpy.nanstd(column_shifts, axis=0) + median_shift = bn.nanmedian(column_shifts, axis=0) + std_shift = bn.nanstd(column_shifts, axis=0) if numpy.abs(median_shift) > 0.5: log.warning(f"large thermal shift measured: {','.join(map(str, column_shifts))} pixels for {mjd = }, {expnum = }, {camera = }") image.add_header_comment(f"large thermal shift: {','.join(map(str, column_shifts))} pixels {camera = }") @@ -413,8 +413,8 @@ def _fix_fiber_thermal_shifts(image, trace_cent, trace_width=None, trace_amp=Non for i, (bmin, bmax) in enumerate(blocks_bounds): x = numpy.arange(bmax-bmin) + i*(bmax-bmin) + 10 # y_models = fiber_model._data[bmin:bmax,c-column_width:c+column_width] - y_model = numpy.nanmedian(fiber_model._data[bmin:bmax,c-column_width:c+column_width], axis=1) - y_data = numpy.nanmedian(image._data[bmin:bmax, c-column_width:c+column_width], axis=1) + y_model = bn.nanmedian(fiber_model._data[bmin:bmax,c-column_width:c+column_width], axis=1) + y_data = bn.nanmedian(image._data[bmin:bmax, c-column_width:c+column_width], axis=1) snr = numpy.sqrt(y_data.mean()) y_model = _normalize_peaks(y_model, min_peak_dist=5.0) y_data = _normalize_peaks(y_data, min_peak_dist=5.0) @@ -1659,7 +1659,7 @@ def find_peaks_auto( ax.plot(pixels, peaks, "o", color="tab:red", mew=0, ms=5) ax.plot( centers, - numpy.ones(len(centers)) * numpy.nanmax(peaks) * 0.5, + numpy.ones(len(centers)) * bn.nanmax(peaks) * 0.5, "x", mew=1, ms=7, @@ -2015,8 +2015,8 @@ def findPeaksMaster2_drp( # find location of peaks (local maxima) either above a fixed threshold or to reach a fixed number of peaks peaks_good = [] - if numpy.nanmax(cut._data) < threshold: - threshold = numpy.nanmax(cut._data) * 0.8 + if bn.nanmax(cut._data) < threshold: + threshold = bn.nanmax(cut._data) * 0.8 while len(peaks_good) != numpy.sum(select_good): (peaks_good, temp, peaks_flux) = cut.findPeaks(threshold=threshold, npeaks=0) if peaks_good[0] < border: @@ -2184,7 +2184,7 @@ def trace_peaks( profile = img.getSlice(ref_column, axis="y")._data if correct_ref: ypix = numpy.arange(profile.size) - guess_heights = numpy.ones_like(positions) * numpy.nanmax(profile) + guess_heights = numpy.ones_like(positions) * bn.nanmax(profile) ref_profile = _spec_from_lines(positions, sigma=1.2, wavelength=ypix, heights=guess_heights) log.info(f"correcting guess positions for column {ref_column}") cc, bhat, mhat = _cross_match( @@ -2654,10 +2654,10 @@ def subtract_straylight( cbar.set_label(f"Counts ({unit})", fontsize="small", color="tab:red") colors_x = plt.cm.coolwarm(numpy.linspace(0, 1, img_median._data.shape[0])) colors_y = plt.cm.coolwarm(numpy.linspace(0, 1, img_median._data.shape[1])) - ax_strayx.fill_between(x_pixels, numpy.nanmedian(img._error, axis=0), lw=0, fc="0.8") + ax_strayx.fill_between(x_pixels, bn.nanmedian(img._error, axis=0), lw=0, fc="0.8") for iy in y_pixels: ax_strayx.plot(x_pixels, img_stray._data[iy], ",", color=colors_x[iy], alpha=0.2) - ax_strayy.fill_betweenx(y_pixels, 0, numpy.nanmedian(img._error, axis=1), lw=0, fc="0.8") + ax_strayy.fill_betweenx(y_pixels, 0, bn.nanmedian(img._error, axis=1), lw=0, fc="0.8") for ix in x_pixels: ax_strayy.plot(img_stray._data[:, ix], y_pixels, ",", color=colors_y[ix], alpha=0.2) save_fig(fig, product_path=out_image, to_display=display_plots, figure_path="qa", label="straylight_model") @@ -2942,7 +2942,7 @@ def offsetTrace_drp( offsets = [] else: offsets.append( - numpy.nanmedian( + bn.nanmedian( numpy.array(log_lines[i + 2].split()[1:]).astype("float32") ) ) @@ -3004,11 +3004,11 @@ def offsetTrace_drp( for j in range(len(out[0])): string_x += " %.3f" % (out[1][j]) string_y += " %.3f" % (out[0][j]) - string_pix += " %.3f" % (numpy.nanmedian(block_line_pos[j])) + string_pix += " %.3f" % (bn.nanmedian(block_line_pos[j])) log.write(string_x + "\n") log.write(string_pix + "\n") log.write(string_y + "\n") - off_trace_median = numpy.nanmedian(numpy.array(off_trace_all)) + off_trace_median = bn.nanmedian(numpy.array(off_trace_all)) off_trace_rms = numpy.std(numpy.array(off_trace_all)) off_trace_rms = "%.4f" % off_trace_rms if numpy.isfinite(off_trace_rms) else "NAN" img.setHdrValue( @@ -3150,12 +3150,12 @@ def offsetTrace2_drp( for j in range(len(out[0])): string_x += " %.3f" % (out[1][j]) string_y += " %.3f" % (out[0][j] * -1) - string_pix += " %.3f" % (numpy.nanmedian(block_line_pos[j])) + string_pix += " %.3f" % (bn.nanmedian(block_line_pos[j])) log.write(string_x + "\n") log.write(string_pix + "\n") log.write(string_y + "\n") - off_trace_median = numpy.nanmedian(numpy.array(off_trace_all)) + off_trace_median = bn.nanmedian(numpy.array(off_trace_all)) off_trace_rms = numpy.std(numpy.array(off_trace_all)) img.setHdrValue( "HIERARCH PIPE FLEX YOFF", @@ -3366,7 +3366,7 @@ def extract_spectra( # propagate thermal shift to slitmap channel = img._header['CCD'][0] slitmap[f"ypix_{channel}"] = slitmap[f"ypix_{channel}"].astype("float64") - slitmap[f"ypix_{channel}"][select_spec] += numpy.nanmedian(shifts, axis=0) + slitmap[f"ypix_{channel}"][select_spec] += bn.nanmedian(shifts, axis=0) if error is not None: error[mask] = replace_error @@ -3385,19 +3385,19 @@ def extract_spectra( rss.setHdrValue("DISPAXIS", 1) rss.setHdrValue( "HIERARCH FIBER CENT MIN", - numpy.nanmin(trace_mask._data[rss._good_fibers]), + bn.nanmin(trace_mask._data[rss._good_fibers]), ) rss.setHdrValue( "HIERARCH FIBER CENT MAX", - numpy.nanmax(trace_mask._data[rss._good_fibers]), + bn.nanmax(trace_mask._data[rss._good_fibers]), ) rss.setHdrValue( "HIERARCH FIBER CENT AVG", - numpy.nanmean(trace_mask._data[rss._good_fibers]) if data.size != 0 else 0, + bn.nanmean(trace_mask._data[rss._good_fibers]) if data.size != 0 else 0, ) rss.setHdrValue( "HIERARCH FIBER CENT MED", - numpy.nanmedian(trace_mask._data[rss._good_fibers]) + bn.nanmedian(trace_mask._data[rss._good_fibers]) if data.size != 0 else 0, ) @@ -3408,19 +3408,19 @@ def extract_spectra( if method == "optimal": rss.setHdrValue( "HIERARCH FIBER WIDTH MIN", - numpy.nanmin(trace_fwhm._data[rss._good_fibers]), + bn.nanmin(trace_fwhm._data[rss._good_fibers]), ) rss.setHdrValue( "HIERARCH FIBER WIDTH MAX", - numpy.nanmax(trace_fwhm._data[rss._good_fibers]), + bn.nanmax(trace_fwhm._data[rss._good_fibers]), ) rss.setHdrValue( "HIERARCH FIBER WIDTH AVG", - numpy.nanmean(trace_fwhm._data[rss._good_fibers]) if data.size != 0 else 0, + bn.nanmean(trace_fwhm._data[rss._good_fibers]) if data.size != 0 else 0, ) rss.setHdrValue( "HIERARCH FIBER WIDTH MED", - numpy.nanmedian(trace_fwhm._data[rss._good_fibers]) + bn.nanmedian(trace_fwhm._data[rss._good_fibers]) if data.size != 0 else 0, ) @@ -3910,8 +3910,8 @@ def preproc_raw_frame( os_models.append(os_model) # compute overscan stats - os_bias_med[i] = numpy.nanmedian(os_quad._data, axis=None) - os_bias_std[i] = numpy.nanmedian(numpy.nanstd(os_quad._data, axis=1), axis=None) + os_bias_med[i] = bn.nanmedian(os_quad._data, axis=None) + os_bias_std[i] = bn.nanmedian(bn.nanstd(os_quad._data, axis=1), axis=None) log.info( f"median and standard deviation in OS quadrant {i+1}: " f"{os_bias_med[i]:.2f} +/- {os_bias_std[i]:.2f} (ADU)" @@ -4048,8 +4048,8 @@ def preproc_raw_frame( axis=0, nstrip=1, ax=axs[i], - mu_stat=numpy.nanmedian, - sg_stat=lambda x, axis: numpy.nanmedian(numpy.std(x, axis=axis)), + mu_stat=bn.nanmedian, + sg_stat=lambda x, axis: bn.nanmedian(numpy.std(x, axis=axis)), labels=True, ) os_x, os_y = _parse_ccd_section(list(os_sections)[0]) @@ -4083,8 +4083,8 @@ def preproc_raw_frame( axis=1, nstrip=1, ax=axs[i], - mu_stat=numpy.nanmedian, - sg_stat=lambda x, axis: numpy.nanmedian(numpy.std(x, axis=axis)), + mu_stat=bn.nanmedian, + sg_stat=lambda x, axis: bn.nanmedian(numpy.std(x, axis=axis)), show_individuals=True, labels=True, ) @@ -4112,8 +4112,8 @@ def preproc_raw_frame( axis=1, nstrip=1, ax=axs[i], - mu_stat=numpy.nanmedian, - sg_stat=lambda x, axis: numpy.nanmedian(numpy.std(x, axis=axis)), + mu_stat=bn.nanmedian, + sg_stat=lambda x, axis: bn.nanmedian(numpy.std(x, axis=axis)), labels=True, ) axs[i].step( @@ -4121,7 +4121,7 @@ def preproc_raw_frame( ) axs[i].step(numpy.arange(os_profiles[i].size), os_models[i], color="k", lw=1) axs[i].axhline( - numpy.nanmedian(os_quad._data.flatten()) + rdnoise[i], + bn.nanmedian(os_quad._data.flatten()) + rdnoise[i], ls="--", color="tab:purple", lw=1, @@ -4436,7 +4436,7 @@ def detrend_frame( quad.computePoissonError(rdnoise) bcorr_img.setSection(section=quad_sec, subimg=quad, inplace=True) - log.info(f"median error in quadrant {i+1}: {numpy.nanmedian(quad._error):.2f} (e-)") + log.info(f"median error in quadrant {i+1}: {bn.nanmedian(quad._error):.2f} (e-)") bcorr_img.setHdrValue("BUNIT", "electron", "physical units of the image") else: @@ -4594,7 +4594,7 @@ def create_master_frame(in_images: List[str], out_image: str, batch_size: int = master_img = combineImages(org_imgs, method="median", normalize=False) elif master_type == "pixflat": master_img = combineImages( - [img / numpy.nanmedian(img._data) for img in org_imgs], + [img / bn.nanmedian(img._data) for img in org_imgs], method="median", normalize=True, normalize_percentile=75, diff --git a/python/lvmdrp/functions/plotMethod.py b/python/lvmdrp/functions/plotMethod.py index ecf683ad..1d5813a2 100644 --- a/python/lvmdrp/functions/plotMethod.py +++ b/python/lvmdrp/functions/plotMethod.py @@ -141,8 +141,8 @@ def flexurePatternTarget_drp(trace_log, disp_log, figure, object, scale="200"): ax.arrow( 500, 1000, - numpy.median(obj_disp_offset[select_disp][0]) * scale, - numpy.median(obj_trace_offset[select_obj][0]) * scale, + bn.median(obj_disp_offset[select_disp][0]) * scale, + bn.median(obj_trace_offset[select_obj][0]) * scale, width=5, head_width=20, head_length=10, diff --git a/python/lvmdrp/functions/rssMethod.py b/python/lvmdrp/functions/rssMethod.py index 0c139574..22144630 100644 --- a/python/lvmdrp/functions/rssMethod.py +++ b/python/lvmdrp/functions/rssMethod.py @@ -97,10 +97,10 @@ def _illumination_correction(fiberflat, apply_correction=True): data[(fiberflat._mask)|(data <= 0)] = numpy.nan # compute median factors - sci_factor = numpy.nanmedian(data[sci_fibers, 1000:3000]) - skw_factor = numpy.nanmedian(data[skw_fibers, 1000:3000]) - ske_factor = numpy.nanmedian(data[ske_fibers, 1000:3000]) - std_factor = numpy.nanmedian(data[std_fibers, 1000:3000]) + sci_factor = bn.nanmedian(data[sci_fibers, 1000:3000]) + skw_factor = bn.nanmedian(data[skw_fibers, 1000:3000]) + ske_factor = bn.nanmedian(data[ske_fibers, 1000:3000]) + std_factor = bn.nanmedian(data[std_fibers, 1000:3000]) norm = numpy.mean([sci_factor, skw_factor, ske_factor, std_factor]) sci_factor /= norm skw_factor /= norm @@ -394,7 +394,7 @@ def determine_wavelength_solution(in_arcs: List[str]|str, out_wave: str, out_lsf if negative: log.info("flipping arc along flux direction") - arc = -1 * arc + numpy.nanmedian(arc._data) + arc = -1 * arc + bn.nanmedian(arc._data) # setup storage array wave_coeffs = numpy.zeros((arc._fibers, numpy.abs(poly_disp) + 1)) @@ -502,12 +502,12 @@ def determine_wavelength_solution(in_arcs: List[str]|str, out_wave: str, out_lsf wave_coeffs[i, :] = wave_poly.convert().coef wave_sol[i, :] = wave_poly(arc._pixels) - wave_rms[i] = numpy.nanstd(wave_poly(cent_wave[i, good_lines]) - ref_lines[good_lines]) + wave_rms[i] = bn.nanstd(wave_poly(cent_wave[i, good_lines]) - ref_lines[good_lines]) log.info( "finished wavelength fitting with median " - f"RMS = {numpy.nanmedian(wave_rms):g} Angstrom " - f"({numpy.nanmedian(wave_rms[:,None]/numpy.diff(wave_sol, axis=1)):g} pix)" + f"RMS = {bn.nanmedian(wave_rms):g} Angstrom " + f"({bn.nanmedian(wave_rms[:,None]/numpy.diff(wave_sol, axis=1)):g} pix)" ) # Estimate the spectral resolution pattern @@ -541,12 +541,12 @@ def determine_wavelength_solution(in_arcs: List[str]|str, out_wave: str, out_lsf lsf_coeffs[i, :] = fwhm_poly.convert().coef lsf_sol[i, :] = fwhm_poly(arc._pixels) - lsf_rms[i] = numpy.nanstd(fwhm_wave - fwhm_poly(cent_wave[i, good_lines])) + lsf_rms[i] = bn.nanstd(fwhm_wave - fwhm_poly(cent_wave[i, good_lines])) log.info( "finished LSF fitting with median " - f"RMS = {numpy.nanmedian(lsf_rms):g} Angstrom " - f"({numpy.nanmedian(lsf_rms[:,None]/numpy.gradient(wave_sol, axis=1)):g} pix)" + f"RMS = {bn.nanmedian(lsf_rms):g} Angstrom " + f"({bn.nanmedian(lsf_rms[:,None]/numpy.gradient(wave_sol, axis=1)):g} pix)" ) # create plot of reference spectrum and wavelength fitting residuals @@ -607,7 +607,7 @@ def determine_wavelength_solution(in_arcs: List[str]|str, out_wave: str, out_lsf ) arc.setHdrValue( "HIERARCH PIPE DISP RMS MEDIAN", - "%.4f" % (numpy.median(wave_rms[good_fibers])), + "%.4f" % (bn.median(wave_rms[good_fibers])), "Median RMS of disp sol", ) arc.setHdrValue( @@ -627,7 +627,7 @@ def determine_wavelength_solution(in_arcs: List[str]|str, out_wave: str, out_lsf ) arc.setHdrValue( "HIERARCH PIPE DISP RMS MEDIAN", - "%.4f" % (numpy.median(lsf_rms[good_fibers])), + "%.4f" % (bn.median(lsf_rms[good_fibers])), "Median RMS of disp sol", ) arc.setHdrValue( @@ -713,7 +713,7 @@ def shift_wave_skylines(in_frame: str, out_frame: str, dwave: float = 8.0, skyli continue offsets[:, ifiber] = sky_wave - skylines - fiber_offset[ifiber] = numpy.nanmedian(offsets[:,ifiber], axis=0) + fiber_offset[ifiber] = bn.nanmedian(offsets[:,ifiber], axis=0) # split per spectrographs specoffset = numpy.asarray(numpy.split(offsets, 3, axis=1)) @@ -912,8 +912,8 @@ def checkPixTable_drp( "%.3f %.3f %.3f %.3f \n" % ( centres[j], - numpy.median(fit_wave[good_fiber, j]), - numpy.median(fit_wave[good_fiber, j]) - centres[j], + bn.median(fit_wave[good_fiber, j]), + bn.median(fit_wave[good_fiber, j]) - centres[j], numpy.std(fit_wave[good_fiber, j]), ) ) @@ -923,7 +923,7 @@ def checkPixTable_drp( for i in range(len(blocks)): if numpy.sum(blocks_good[i]) > 0: log.write( - " %.3f" % numpy.median(offset_pix[blocks[i][blocks_good[i]], j]) + " %.3f" % bn.median(offset_pix[blocks[i][blocks_good[i]], j]) ) else: log.write(" 0.0") @@ -933,7 +933,7 @@ def checkPixTable_drp( log.write( " %.3f" % ( - numpy.median(fit_wave[blocks[i][blocks_good[i]], j]) + bn.median(fit_wave[blocks[i][blocks_good[i]], j]) - centres[j] ) ) @@ -941,7 +941,7 @@ def checkPixTable_drp( log.write(" 0.0") log.write("\n") - off_disp_median = numpy.median(offset_pix[good_fiber, :]) + off_disp_median = bn.median(offset_pix[good_fiber, :]) off_disp_rms = numpy.std(offset_pix[good_fiber, :]) off_disp_median = ( float("%.4f" % off_disp_median) @@ -1060,11 +1060,11 @@ def correctPixTable_drp( for i in range(rss._fibers): spec = rss[i] if smooth_poly_disp == "": - off = numpy.median(offsets.flatten()) + off = bn.median(offsets.flatten()) else: smooth_poly_disp = int(smooth_poly_disp) if smooth_poly_disp == "": - off = numpy.median(offsets[i]) + off = bn.median(offsets[i]) else: off = Spectrum1D(wave=ref_wave, data=offsets[:, i]) off.smoothPoly(smooth_poly_disp, ref_base=spec._wave) @@ -1518,7 +1518,7 @@ def correctTraceMask_drp(trace_in, trace_out, logfile, ref_file, poly_smooth="") trace = TraceMask.from_file(trace_in) if poly_smooth == "": - trace = trace + (numpy.median(offsets.flatten()) * -1) + trace = trace + (bn.median(offsets.flatten()) * -1) else: split_trace = trace.split(offsets.shape[1], axis="y") offset_trace = TraceMask() @@ -1788,7 +1788,7 @@ def matchFluxRSS_drp( # load subimages from disc and append them to a list rss = loadRSS(list_rss[i]) specs.append(rss.createAperSpec(center_x, center_y, arc_radius)) - fluxes.append(numpy.median(specs[i]._data)) + fluxes.append(bn.median(specs[i]._data)) order = numpy.argsort(fluxes) # print fluxes, order diff --git a/python/lvmdrp/functions/run_calseq.py b/python/lvmdrp/functions/run_calseq.py index bb3f694e..b72483ba 100644 --- a/python/lvmdrp/functions/run_calseq.py +++ b/python/lvmdrp/functions/run_calseq.py @@ -1401,7 +1401,7 @@ def create_dome_fiberflats(mjd, expnums_ldls, expnums_qrtz, use_longterm_cals=Tr # normalize by median fiber fflat = RSS(data=mamp._data, error=np.sqrt(mamp._data), mask=xflat._mask, wave_trace=mwave, lsf_trace=mlsf, header=xflat._header) fflat = fflat.rectify_wave(method="linear", wave_range=SPEC_CHANNELS[channel], wave_disp=0.5) - median_fiber = np.nanmedian(fflat._data, axis=0) + median_fiber = bn.nanmedian(fflat._data, axis=0) fflat._data = fflat._data / median_fiber fflat.set_wave_trace(mwave) fflat.set_lsf_trace(mlsf) diff --git a/python/lvmdrp/functions/run_twilights.py b/python/lvmdrp/functions/run_twilights.py index 0bfd632f..dffb0498 100644 --- a/python/lvmdrp/functions/run_twilights.py +++ b/python/lvmdrp/functions/run_twilights.py @@ -16,6 +16,7 @@ from mpl_toolkits.axes_grid1.inset_locator import inset_axes from matplotlib.gridspec import GridSpec +import bottleneck as bn from lvmdrp import log from lvmdrp.core.constants import SPEC_CHANNELS from lvmdrp.core.tracemask import TraceMask @@ -92,7 +93,7 @@ def remove_field_gradient(in_hflat, out_gflat, wrange, deg=1, display_plots=Fals x, y, flux, grad_model, grad_model_e, grad_model_w, grad_model_s = fflat.fit_field_gradient(wrange, poly_deg=deg) flux_c = flux / grad_model - flux_med = np.nanmedian(flux) + flux_med = bn.nanmedian(flux) flux_std = np.nanstd(flux) fig, axs = create_subplots(to_display=display_plots, nrows=1, ncols=3, figsize=(15,6), sharex=True, sharey=True, layout="constrained") @@ -345,8 +346,8 @@ def fit_fiberflat(in_twilight: str, out_flat: str, out_twilight: str, axs[1].set_xlim(*axs[0].get_xlim()) ypixels = np.split(np.arange(flat_twilight._fibers) + 1, 3) - median = np.split(np.nanmedian(flat_twilight._data, axis=1), 3) - mu = np.nanmedian(flat_twilight._data) + median = np.split(bn.nanmedian(flat_twilight._data, axis=1), 3) + mu = bn.nanmedian(flat_twilight._data) median = (median / mu - 1) * 100 axs[2].axhspan(-0.5, 0.5, lw=0, alpha=0.3, color="0.7") axs[2].axhline(0, ls="--", lw=1, color="0.7") @@ -474,7 +475,7 @@ def combine_twilight_sequence(in_fiberflats: List[str], out_fiberflat: str, # get exposed standard fiber ID fiber_id = fflat._header.get("CALIBFIB") if fiber_id is None: - snrs = np.nanmedian(fflat._data / fflat._error, axis=1) + snrs = bn.nanmedian(fflat._data / fflat._error, axis=1) select_nonexposed = snrs < 50 #plt.figure(figsize=(15,5)) #plt.plot(snrs[select_allstd]) From 2afb42b5f1a9dc2f7d7268bfc0e1b845272bca39 Mon Sep 17 00:00:00 2001 From: Niv Drory Date: Tue, 10 Sep 2024 10:05:33 -0500 Subject: [PATCH 2/5] use bottleneck instead of numpy median/nanmedian --- python/lvmdrp/functions/specialMethod.py | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/python/lvmdrp/functions/specialMethod.py b/python/lvmdrp/functions/specialMethod.py index e19bfd1d..a1985702 100644 --- a/python/lvmdrp/functions/specialMethod.py +++ b/python/lvmdrp/functions/specialMethod.py @@ -1,6 +1,7 @@ import numpy import pylab +import bottleneck as bn from lvmdrp.core.fiberrows import FiberRows from lvmdrp.core.fit_profile import Exponential_constant from lvmdrp.core.header import Header @@ -108,14 +109,14 @@ def extinctCAVEX_drp( dates == date_in, numpy.logical_and(UThours >= 0.0, UThours <= 14.0) ) if numpy.sum(select_time) > 0: - out_av = numpy.median(av[select_time]) + out_av = bn.median(av[select_time]) else: select_time = numpy.logical_or( numpy.logical_and(dates == date_in, UThours > 14.0), numpy.logical_and(dates == date_in + 1, UThours < 14.0), ) if numpy.sum(select_time) > 0: - out_av = numpy.median(av[select_time]) + out_av = bn.median(av[select_time]) else: out_av = missing_extinct else: @@ -131,7 +132,7 @@ def extinctCAVEX_drp( numpy.logical_and(dates == date_max, UThours < 14.0), ) if numpy.sum(select_time) > 0: - out_av = numpy.median(av[select_time]) + out_av = bn.median(av[select_time]) else: out_av = missing_extinct @@ -383,14 +384,14 @@ def checkWavelengthRSS_drp( % ( line_name[i], line_list[i], - numpy.median(cent_wave[mask_line, i]) * cdelt + crval, - numpy.mean(cent_wave[mask_line, i]) * cdelt + crval, - numpy.std(cent_wave[mask_line, i]) * cdelt, - numpy.median(flux[mask_line, i]), - numpy.std(flux[mask_line, i]), - numpy.median(fwhm[mask_line, i]) * cdelt, - numpy.mean(fwhm[mask_line, i]) * cdelt, - numpy.std(fwhm[mask_line, i]) * cdelt, + bn.median(cent_wave[mask_line, i]) * cdelt + crval, + bn.mean(cent_wave[mask_line, i]) * cdelt + crval, + bn.std(cent_wave[mask_line, i]) * cdelt, + bn.median(flux[mask_line, i]), + bn.std(flux[mask_line, i]), + bn.median(fwhm[mask_line, i]) * cdelt, + bn.mean(fwhm[mask_line, i]) * cdelt, + bn.std(fwhm[mask_line, i]) * cdelt, ) ) From 08466c9e0e4af0fdb841f1dc2eda44d2cb4aba82 Mon Sep 17 00:00:00 2001 From: Niv Drory Date: Tue, 10 Sep 2024 13:47:03 -0500 Subject: [PATCH 3/5] add missing bottleneck import --- python/lvmdrp/functions/run_calseq.py | 1 + 1 file changed, 1 insertion(+) diff --git a/python/lvmdrp/functions/run_calseq.py b/python/lvmdrp/functions/run_calseq.py index b72483ba..8036dbb1 100644 --- a/python/lvmdrp/functions/run_calseq.py +++ b/python/lvmdrp/functions/run_calseq.py @@ -27,6 +27,7 @@ import os import numpy as np +import bottleneck as bn from glob import glob from copy import deepcopy as copy from shutil import copy2, rmtree From 3c6fbe98fb517dc82e1c07d8dc982f0c6d737db7 Mon Sep 17 00:00:00 2001 From: Niv Drory Date: Tue, 10 Sep 2024 13:53:21 -0500 Subject: [PATCH 4/5] add missing bottleneck imports --- python/lvmdrp/functions/cubeMethod.py | 2 +- python/lvmdrp/functions/plotMethod.py | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/python/lvmdrp/functions/cubeMethod.py b/python/lvmdrp/functions/cubeMethod.py index c5516276..ddb77383 100644 --- a/python/lvmdrp/functions/cubeMethod.py +++ b/python/lvmdrp/functions/cubeMethod.py @@ -1,7 +1,7 @@ import sys import numpy - +import bottleneck as bn try: import pylab diff --git a/python/lvmdrp/functions/plotMethod.py b/python/lvmdrp/functions/plotMethod.py index 1d5813a2..1e8cba40 100644 --- a/python/lvmdrp/functions/plotMethod.py +++ b/python/lvmdrp/functions/plotMethod.py @@ -5,6 +5,7 @@ except ImportError: pass import numpy +import bottleneck as bn description = "Provides Methods to make some plots" From f18f1b51196e8a0d93e87bec245f7a2c8a99c325 Mon Sep 17 00:00:00 2001 From: Alfredo Mejia-Narvaez Date: Tue, 10 Sep 2024 16:10:39 -0300 Subject: [PATCH 5/5] adding missing import --- python/lvmdrp/core/cube.py | 1 + 1 file changed, 1 insertion(+) diff --git a/python/lvmdrp/core/cube.py b/python/lvmdrp/core/cube.py index 6260090c..e94c3685 100644 --- a/python/lvmdrp/core/cube.py +++ b/python/lvmdrp/core/cube.py @@ -1,6 +1,7 @@ import numpy from astropy.io import fits as pyfits from scipy import ndimage +import bottleneck as bn from lvmdrp.core.header import Header, combineHdr from lvmdrp.core.image import Image