Merge branch 'master' of github.com:gbrammer/grizli

grizli/fitting.py

@@ -1172,7 +1172,7 @@ def compute_cdf_percentiles(fit, cdf_sigmas=CDF_SIGMAS):
                         
     """
     from scipy.interpolate import Akima1DInterpolator
-    from scipy.integrate import cumtrapz
+    from scipy.integrate import cumulative_trapezoid
     import scipy.stats
 
     if cdf_sigmas is None:
@@ -1190,7 +1190,7 @@ def compute_cdf_percentiles(fit, cdf_sigmas=CDF_SIGMAS):
     pz_fine = np.exp(spl(zfine))
     pz_fine[~ok] = 0
 
-    cdf_fine = cumtrapz(pz_fine, x=zfine)
+    cdf_fine = cumulative_trapezoid(pz_fine, x=zfine)
     cdf_x = np.interp(cdf_y, cdf_fine/cdf_fine[-1], zfine[1:])
 
     return cdf_x, cdf_y
@@ -2513,10 +2513,9 @@ def xfit_redshift(self, prior=None,
             
             
         """
-        from numpy.polynomial.polynomial import polyfit, polyval
+        from numpy.polynomial import Polynomial
         from scipy.stats import t as student_t
         from scipy.special import huber
-        import peakutils
 
         if isinstance(zr, int):
             if zr == 0:
@@ -2657,7 +2656,7 @@ def xfit_redshift(self, prior=None,
 
         if len(zgrid) > 1:
             chi2_rev[chi2_rev < 0] = 0
-            indexes = peakutils.indexes(chi2_rev, thres=0.4, min_dist=8)
+            indexes = utils.find_peaks(chi2_rev,threshold=0.4,min_dist=9)
             num_peaks = len(indexes)
             so = np.argsort(chi2_rev[indexes])
             indexes = indexes[so[::-1]]
@@ -2676,9 +2675,9 @@ def xfit_redshift(self, prior=None,
             zgrid_zoom = []
             for ix in indexes[:max_peaks]:
                 if (ix > 0) & (ix < len(chi2)-1):
-                    c = polyfit(zgrid[ix-1:ix+2], chi2[ix-1:ix+2], 2)
-                    zi = -c[1]/(2*c[0])
-                    chi_i = polyval(c, zi)
+                    p = Polynomial.fit(zgrid[ix-1:ix+2], chi2[ix-1:ix+2], deg=2)
+                    zi = p.deriv().roots()[0]
+                    chi_i = p(zi)
 
                     # Just use grid value
                     zi = zgrid[ix]
@@ -2857,9 +2856,10 @@ def _parse_zfit_output(self, fit, risk_gamma=0.15, prior=None):
         Adds metadata and `pdf` and `risk` columns to `fit` table
         
         """
+        from numpy.polynomial import Polynomial
         import scipy.interpolate
         from scipy.interpolate import Akima1DInterpolator
-        from scipy.integrate import cumtrapz
+        from scipy.integrate import cumulative_trapezoid
 
         # Normalize to min(chi2)/DoF = 1.
         scl_nu = fit['chi2'].min()/self.DoF
@@ -2896,7 +2896,7 @@ def _parse_zfit_output(self, fit, risk_gamma=0.15, prior=None):
             # Compute CDF and probability intervals
             #dz = np.gradient(zfine[ok])
             #cdf = np.cumsum(np.exp(spl(zfine[ok]))*dz/norm)
-            cdf = cumtrapz(pz_fine, x=zfine)
+            cdf = cumulative_trapezoid(pz_fine, x=zfine)
             percentiles = np.array([2.5, 16, 50, 84, 97.5])/100.
             pz_percentiles = np.interp(percentiles, cdf/cdf[-1], zfine[1:])
 
@@ -2912,8 +2912,8 @@ def _parse_zfit_output(self, fit, risk_gamma=0.15, prior=None):
             # Fit a parabola around min(risk)
             zi = np.argmin(risk)
             if (zi < len(risk)-1) & (zi > 0):
-                c = np.polyfit(fit['zgrid'][zi-1:zi+2], risk[zi-1:zi+2], 2)
-                z_risk = -c[1]/(2*c[0])
+                p = Polynomial.fit(fit['zgrid'][zi-1:zi+2], risk[zi-1:zi+2],deg=2)
+                z_risk = p.deriv().roots()[0]
             else:
                 z_risk = fit['zgrid'][zi]
 
@@ -2924,8 +2924,8 @@ def _parse_zfit_output(self, fit, risk_gamma=0.15, prior=None):
             # MAP, maximum p(z) from parabola fit around tabulated maximum
             zi = np.argmax(pdf)
             if (zi < len(pdf)-1) & (zi > 0):
-                c = np.polyfit(fit['zgrid'][zi-1:zi+2], pdf[zi-1:zi+2], 2)
-                z_map = -c[1]/(2*c[0])
+                p = Polynomial.fit(fit['zgrid'][zi-1:zi+2], pdf[zi-1:zi+2],deg=2)
+                z_map = p.deriv().roots()[0]
             else:
                 z_map = fit['zgrid'][zi]
         else:
@@ -3449,8 +3449,7 @@ def compute_scale_array(pscale, wave):
         pscale : array-like
             Coefficients of the linear model normalized by factors of 10 per 
             order, i.e, ``pscale = [10]`` is a constant unit scaling.  Note 
-            that parameter order is reverse that expected by 
-            `numpy.polyval`.
+            that parameter is what is expected by `numpy.polynomial.Polynomial`.
         
         wave : array-like
             Wavelength grid in Angstroms.  Scaling is normalized to 
@@ -3464,12 +3463,12 @@ def compute_scale_array(pscale, wave):
             >>> pscale = [10]
             >>> N = len(pscale)
             >>> rescale = 10**(np.arange(N)+1)
-            >>> wscale = np.polyval((pscale/rescale)[::-1], (wave-1.e4)/1000.)
+            >>> wscale = np.polynomial.Polynomial(pscale/rescale)((wave-1.e4)/1000.)
 
         """
         N = len(pscale)
         rescale = 10**(np.arange(N)+1)
-        wscale = np.polyval((pscale/rescale)[::-1], (wave-1.e4)/1000.)
+        wscale = np.polynomial.Polynomial(pscale/rescale)((wave-1.e4)/1000.)
         return wscale
 
 
@@ -3480,7 +3479,7 @@ def objfun_scale(pscale, AxT, data, self, retval):
         spectra
         """
         import scipy.optimize
-        from numpy.polynomial.polynomial import polyval
+        from numpy.polynomial import Polynomial
 
         scale = self.compute_scale_array(pscale, self.wavef[self.fit_mask])
         scale[-self.Nphot:] = 1.
@@ -4645,7 +4644,7 @@ def get_flat_model(self, spectrum_1d, id=None, apply_mask=True, is_cgs=True):
         """
         mfull = []
         for ib, beam in enumerate(self.beams):
-            if spectrum_1d is -1:
+            if spectrum_1d == -1:
                 model_i = beam.model*1
             else:
                 model_i = beam.compute_model(id=id, spectrum_1d=spectrum_1d,

grizli/grismconf.py

@@ -1741,9 +1741,9 @@ def _eval_model(self, model, x0, y0, t, get_coeffs=False):
                     _c.append(m(x0, y0))
 
             if get_coeffs:
-                value = _c[::-1]
+                value = _c
             else:
-                value = np.polyval(_c[::-1], t)
+                value = np.polynomial.Polynomial(_c)(t)
 
         return value
 
@@ -1752,8 +1752,8 @@ def _root_inverse_model(self, model, x0, y0, dx, **kwargs):
         """Inverse values interpolated from the forward model using polynomial roots"""
         coeffs = self._eval_model(model, x0, y0, 0, get_coeffs=True)
         if hasattr(coeffs, '__len__'):
-            coeffs[-1] -= dx
-            value = np.roots(coeffs)[-1]
+            coeffs[0] -= dx
+            value = np.polynomial.Polynomial(coeffs).roots()[-1]
         else:
             value = coeffs
 

grizli/jwst_utils.py

@@ -1024,6 +1024,7 @@ def exposure_oneoverf_correction(file, axis=None, thresholds=[5,4,3], erode_mask
         The row- or column-average correction array
     """
     import numpy as np
+    from numpy.polynomial import Chebyshev
     import scipy.ndimage as nd
     import matplotlib.pyplot as plt
 
@@ -1150,8 +1151,8 @@ def exposure_oneoverf_correction(file, axis=None, thresholds=[5,4,3], erode_mask
         deg = nx//deg_pix
 
         for _iter in range(3):
-            coeffs = np.polynomial.chebyshev.chebfit(xarr[ok], med[ok], deg)
-            cheb = np.polynomial.chebyshev.chebval(xarr, coeffs)
+            cfit = Chebyshev.fit(xarr[ok], med[ok], deg=deg)
+            cheb = cfit(xarr)
             ok = np.isfinite(med) & (np.abs(med-cheb) < 0.05)
 
         if make_plot:
@@ -1408,7 +1409,7 @@ def initialize_jwst_image(filename, verbose=True, max_dq_bit=14, orig_keys=ORIG_
             try:
                 _ = exposure_oneoverf_correction(filename, in_place=True, 
                                              **oneoverf_kwargs)
-            except TypeError:
+            except ValueError:
                 # Should only fail for test data
                 utils.log_exception(utils.LOGFILE, traceback)
                 msg = f'exposure_oneoverf_correction: failed for {filename}'
@@ -1422,7 +1423,7 @@ def initialize_jwst_image(filename, verbose=True, max_dq_bit=14, orig_keys=ORIG_
                         _ = exposure_oneoverf_correction(filename, in_place=True, 
                                                          axis=-1,
                                                          **oneoverf_kwargs)
-                    except TypeError:
+                    except ValueError:
                         # Should only fail for test data
                         utils.log_exception(utils.LOGFILE, traceback)
                         msg = f'exposure_oneoverf_correction: axis=-1 failed for {filename}'

grizli/model.py

@@ -405,7 +405,7 @@ def process_config(self):
         self.NX = len(self.dx)
         self.sh_beam = (self.sh[0], self.sh[1]+self.NX)
 
-        self.modelf = np.zeros(np.product(self.sh_beam), dtype=np.float32)
+        self.modelf = np.zeros(np.prod(self.sh_beam), dtype=np.float32)
         self.model = self.modelf.reshape(self.sh_beam)
         self.idx = np.arange(self.modelf.size,
                              dtype=np.int64).reshape(self.sh_beam)
@@ -3260,7 +3260,7 @@ def compute_model_orders(self,
             # if psf_params is not None:
             #     skip_init_psf = False
             #     if hasattr(beam, 'psf_params'):
-            #         skip_init_psf |= np.product(np.isclose(beam.psf_params, psf_params)) > 0
+            #         skip_init_psf |= np.prod(np.isclose(beam.psf_params, psf_params)) > 0
             #
             #     if not skip_init_psf:
             #         beam.x_init_epsf(flat_sensitivity=False, psf_params=psf_params, psf_filter=psf_filter, yoff=0.06)
@@ -4758,10 +4758,10 @@ def full_2d_wcs(self, data=None):
         #wcs_header = grizli.utils.to_header(self.grism.wcs)
 
         x = np.arange(len(self.beam.lam_beam))
-        c = np.polyfit(x, self.beam.lam_beam/1.e4, 2)
-        #c = np.polyfit((self.beam.lam_beam-self.beam.lam_beam[0])/1.e4, x/h['CD1_1'], 2)
+        c = np.polynomial.Polynomial.fit(x,self.beam.lam_beam/1.e4,deg=2).convert().coef[::-1]
+        #c = np.polynomial.Polynomial.fit((self.beam.lam_beam-self.beam.lam_beam[0])/1.e4, x/h['CD1_1'],deg=2).convert().coef[::-1]
 
-        ct = np.polyfit(x, self.beam.ytrace_beam, 2)
+        ct = np.polynomial.Polynomial.fit(x,self.beam.ytrace_beam,deg=2).convert().coef[::-1]
 
         h['A_ORDER'] = 2
         h['B_ORDER'] = 2
@@ -4874,7 +4874,7 @@ def get_dispersion_PA(self, decimals=0, local=False):
         if decimals is not None:
             dispersion_PA = np.round(dispersion_PA, decimals=decimals)
 
-        return float(dispersion_PA)
+        return dispersion_PA[0]
 
 
     def init_epsf(self, center=None, tol=1.e-3, yoff=0., skip=1., flat_sensitivity=False, psf_params=None, N=4, get_extended=False, only_centering=True):
@@ -5255,9 +5255,8 @@ def init_poly_coeffs(self, poly_order=1, fit_background=True):
     #             print(utils.NO_NEWLINE + '{0:.4f} {1:9.1f}'.format(zgrid[i], chi2[i]))
     #
     #     # peaks
-    #     import peakutils
     #     chi2nu = (chi2.min()-chi2)/self.DoF
-    #     indexes = peakutils.indexes((chi2nu+0.01)*(chi2nu > -0.004), thres=0.003, min_dist=20)
+    #     indexes = utils.find_peaks((chi2nu+0.01)*(chi2nu > -0.004), threshold=0.003, min_dist=21)
     #     num_peaks = len(indexes)
     #     # plt.plot(zgrid, (chi2-chi2.min())/ self.DoF)
     #     # plt.scatter(zgrid[indexes], (chi2-chi2.min())[indexes]/ self.DoF, color='r')

grizli/multifit.py

@@ -870,7 +870,7 @@ def refine(self, id, mag=-99, poly_order=3, size=30, ds9=None, verbose=True, max
         # Check where templates inconsistent with broad-band fluxes
         xb = [beam.direct.ref_photplam if beam.direct['REF'] is not None else beam.direct.photplam for beam in beams]
         obs_flux = np.array([beam.beam.total_flux for beam in beams])
-        mod_flux = np.polyval(scale_coeffs[::-1], np.array(xb)/1.e4-1)
+        mod_flux = np.polynomial.Polynomial(scale_coeffs)(np.array(xb)/1.e4-1)
         nonz = obs_flux != 0
 
         if (np.abs(mod_flux/obs_flux)[nonz].max() > max_coeff) | ((~np.isfinite(mod_flux/obs_flux)[nonz]).sum() > 0) | (np.min(mod_flux[nonz]) < 0) | ((~np.isfinite(ypoly)).sum() > 0):
@@ -920,7 +920,7 @@ def refine(self, id, mag=-99, poly_order=3, size=30, ds9=None, verbose=True, max
     #     xb = [beam.direct.ref_photplam if beam.direct['REF'] is not None
     #           else beam.direct.photplam for beam in beams]
     #     fb = [beam.beam.total_flux for beam in beams]
-    #     mb = np.polyval(scale_coeffs[::-1], np.array(xb)/1.e4-1)
+    #     mb = np.polynomial.Polynomial(scale_coeffs)(np.array(xb)/1.e4-1)
     #
     #     if (np.abs(mb/fb).max() > max_coeff) | (~np.isfinite(mb/fb).sum() > 0) | (np.min(mb) < 0):
     #         if verbose:
@@ -1860,7 +1860,7 @@ def _parse_beam_arrays(self):
         self.poly_order = None
 
         self.shapes = [beam.model.shape for beam in self.beams]
-        self.Nflat = [np.product(shape) for shape in self.shapes]
+        self.Nflat = [np.prod(shape) for shape in self.shapes]
         self.Ntot = np.sum(self.Nflat)
 
         for b in self.beams:
@@ -2332,7 +2332,7 @@ def eval_poly_spec(self, coeffs_full):
         scale_coeffs = coeffs_full[i0:i0+self.n_poly]
 
         #yspec = [xspec**o*scale_coeffs[o] for o in range(self.poly_order+1)]
-        yfull = np.polyval(scale_coeffs[::-1], xspec - px0)
+        yfull = np.polynomial.Polynomial(scale_coeffs)(xspec - px0)
         return xspec, yfull
 
 
@@ -2751,7 +2751,7 @@ def fit_redshift(self, prior=None, poly_order=1, fwhm=1200,
                      fsps_templates=False):
         """TBD
         """
-        from numpy.polynomial.polynomial import polyfit, polyval
+        from numpy.polynomial import Polynomial
 
         if zr is None:
             zr = [0.65, 1.6]
@@ -2812,16 +2812,15 @@ def fit_redshift(self, prior=None, poly_order=1, fwhm=1200,
         print('First iteration: z_best={0:.4f}\n'.format(zgrid[iz]))
 
         # peaks
-        import peakutils
         # chi2nu = (chi2.min()-chi2)/self.DoF
-        # indexes = peakutils.indexes((chi2nu+delta_chi2_threshold)*(chi2nu > -delta_chi2_threshold), thres=0.3, min_dist=20)
+        # indexes = utils.find_peaks((chi2nu+delta_chi2_threshold)*(chi2nu > - delta_chi2_threshold), threshold=0.3, min_dist=21)
 
         chi2_rev = (chi2_poly - chi2)/self.DoF
         if chi2_poly < (chi2.min() + 9):
             chi2_rev = (chi2.min() + 16 - chi2)/self.DoF
 
         chi2_rev[chi2_rev < 0] = 0
-        indexes = peakutils.indexes(chi2_rev, thres=0.4, min_dist=8)
+        indexes = utils.find_peaks(chi2_rev, threshold=0.4, min_dist=9)
         num_peaks = len(indexes)
 
         if False:
@@ -2834,9 +2833,9 @@ def fit_redshift(self, prior=None, poly_order=1, fwhm=1200,
             zgrid_zoom = []
             for ix in indexes:
                 if (ix > 0) & (ix < len(chi2)-1):
-                    c = polyfit(zgrid[ix-1:ix+2], chi2[ix-1:ix+2], 2)
-                    zi = -c[1]/(2*c[0])
-                    chi_i = polyval(c, zi)
+                    p = Polynomial.fit(zgrid[ix-1:ix+2], chi2[ix-1:ix+2], deg=2)
+                    zi = p.deriv().roots()[0]
+                    chi_i = p(zi)
                     zgrid_zoom.extend(np.arange(zi-2*dz[0],
                                       zi+2*dz[0]+dz[1]/10., dz[1]))
 
@@ -2891,9 +2890,9 @@ def fit_redshift(self, prior=None, poly_order=1, fwhm=1200,
 
         # Fit parabola
         if (ix > 0) & (ix < len(chi2)-1):
-            c = polyfit(zgrid[ix-1:ix+2], chi2[ix-1:ix+2], 2)
-            zbest = -c[1]/(2*c[0])
-            chibest = polyval(c, zbest)
+            p = Polynomial.fit(zgrid[ix-1:ix+2], chi2[ix-1:ix+2], deg=2)
+            zbest = p.deriv().roots()[0]
+            chibest = p(zbest)
 
         out = self.fit_at_z(z=zbest, templates=templates,
                             fitter=fitter, poly_order=poly_order,

grizli/prep.py

@@ -4066,13 +4066,13 @@ def separate_chip_sky(visit, filters=['F200LP','F350LP','F600LP','F390W'], steps
             #deg = 11
 
             for _iter in range(5):
-                cc = np.polynomial.chebyshev.chebfit(xi[msk],
-                                                     row_avg[ext][msk],
-                                                     average_order, 
-                                                     rcond=None,
-                                                     full=False, w=None)
-                                                     
-                row_model[ext] = np.polynomial.chebyshev.chebval(xi, cc)
+                cfit = Chebyshev.fit(xi[msk],
+                                    row_avg[ext][msk],
+                                    deg=average_order,
+                                    rcond=None,
+                                    full=False,
+                                    w=None)
+                row_model[ext] = cfit(xi)            
                 msk = np.isfinite(row_avg[ext]) 
                 msk &= np.abs(row_avg[ext] - row_model[ext]) < 3*utils.nmad(row_avg[ext][msk])
 
@@ -5156,7 +5156,7 @@ def tweak_align(direct_group={}, grism_group={}, max_dist=1., n_min=10, key=' ',
                                  'prep.tweak_align')
 
     from drizzlepac.astrodrizzle import AstroDrizzle
-    from numpy.polynomial.polynomial import polyfit, polyval
+    from numpy.polynomial import Polynomial
 
     if len(direct_group['files']) < 2:
         logstr = '# ! {0}: Only one direct image found, can\'t compute shifts'
@@ -5217,10 +5217,10 @@ def tweak_align(direct_group={}, grism_group={}, max_dist=1., n_min=10, key=' ',
 
         shifts = np.array([shift_dict[k][:2] for k in sorted(shift_dict)])
         t = np.arange(shifts.shape[0])
-        cx = polyfit(t, shifts[:, 0], fit_order)
-        sx = polyval(cx, t)
-        cy = polyfit(t, shifts[:, 1], fit_order)
-        sy = polyval(cy, t)
+        px = Polynomial.fit(t, shifts[:, 0], deg=fit_order)
+        sx = px(t)
+        py = Polynomial.fit(t, shifts[:, 1], deg=fit_order)
+        sy = py(t)
         fit_shift = np.array([sx, sy]).T
 
         for ik, k in enumerate(sorted(shift_dict)):