From 198eb36da6b0bb78ae2a18c210a463bf77c045e9 Mon Sep 17 00:00:00 2001 From: TheSkyentist Date: Thu, 7 Mar 2024 10:09:58 +0100 Subject: [PATCH] Updating to Polynomial --- grizli/fitting.py | 26 +++++++++++++------------- grizli/grismconf.py | 8 ++++---- grizli/model.py | 6 +++--- grizli/multifit.py | 20 ++++++++++---------- grizli/stack.py | 14 +++++++------- 5 files changed, 37 insertions(+), 37 deletions(-) diff --git a/grizli/fitting.py b/grizli/fitting.py index 0394d640..424388a3 100644 --- a/grizli/fitting.py +++ b/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. diff --git a/grizli/grismconf.py b/grizli/grismconf.py index 03870721..ec2ee770 100644 --- a/grizli/grismconf.py +++ b/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 diff --git a/grizli/model.py b/grizli/model.py index 8471f06f..0da82d2a 100644 --- a/grizli/model.py +++ b/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 diff --git a/grizli/multifit.py b/grizli/multifit.py index 83f84706..fbfdcbe0 100644 --- a/grizli/multifit.py +++ b/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, diff --git a/grizli/stack.py b/grizli/stack.py index 7dd84eaf..fcb0e8eb 100755 --- a/grizli/stack.py +++ b/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