Merge pull request #148 from sdss/bottleneck

Consistently use bottleneck

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
@@ -412,7 +413,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

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(

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):

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)

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":
@@ -3484,7 +3484,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)

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,15 +3508,15 @@ 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
                 )
             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:

python/lvmdrp/functions/cubeMethod.py

@@ -1,7 +1,7 @@
 import sys
 
 import numpy
-
+import bottleneck as bn
 
 try:
     import pylab
@@ -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] = (