Updating to Polynomial

grizli/fitting.py

@@ -2513,7 +2513,7 @@ 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
@@ -2676,9 +2676,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,6 +2857,7 @@ 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
@@ -2912,8 +2913,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 +2925,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 +3450,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 +3464,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 +3480,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.

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/model.py

@@ -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

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:
@@ -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]
@@ -2834,9 +2834,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 +2891,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/stack.py

@@ -561,10 +561,10 @@ def scale_AxT(p, Ax, spec_wave, Nphot, Next):
         """
         Scale spectrum templates by polynomial function
         """
-        from numpy.polynomial.polynomial import polyval
+        from np.polynomial import Polynomial
 
         scale = np.ones(Ax.shape[1])
-        scale[:-Nphot] = polyval(p[::-1]/10., (spec_wave-1.e4)/1000.)
+        scale[:-Nphot] = Polynomial(p/10.)((spec_wave-1.e4)/1000.)
         AxT = Ax*scale
         for i in range(Next):
             AxT[i, :] /= scale
@@ -578,10 +578,10 @@ def objective_scale(p, Ax, data, spec_wave, fit_mask, sivarf, Nphot, Next, retur
         spectra
         """
         import scipy.optimize
-        from numpy.polynomial.polynomial import polyval
+        from np.polynomial import Polynomial
 
         scale = np.ones(Ax.shape[1])
-        scale[:-Nphot] = polyval(p[::-1]/10., (spec_wave-1.e4)/1000.)
+        scale[:-Nphot] = Polynomial(p/10.)((spec_wave-1.e4)/1000.)
         AxT = Ax*scale
 
         # Remove scaling from background component
@@ -1422,9 +1422,9 @@ def _build_model(self):
                 # Linear fit to differential G800L dispersion
                 # y = np.diff(beam.wave)/np.diff(beam.wave)[0]
                 # x = beam.wave[1:]
-                # c_i = np.polyfit(x/1.e4, y, 1)
-                c_i = np.array([0.07498747,  0.98928126])
-                scale = np.polyval(c_i, self.wave/1.e4)
+                # c_i = np.polynomial.Polynomial.fit(x/1.e4, y, deg=1)
+                c_i = np.array([0.98928126, 0.07498747])
+                scale = np.polynomial.Polynomial(c_i)(self.wave/1.e4)
                 sens *= scale
 
             self.sens = sens*dlam  # *1.e-17