Skip to content

Commit

Permalink
more joint fit
Browse files Browse the repository at this point in the history
  • Loading branch information
@moustakas
moustakas committed Dec 31, 2022
1 parent 7754c44 commit 5f60414
Show file tree
Hide file tree
Showing 3 changed files with 51 additions and 31 deletions.
39 changes: 23 additions & 16 deletions py/fastspecfit/continuum.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1422,6 +1422,7 @@ def init_spec_output(self, nobj=1):
nssp_coeff = len(self.sspinfo)

out = Table()
out.add_column(Column(name='APCORR', length=nobj, dtype='f4')) # aperture correction
out.add_column(Column(name='CONTINUUM_Z', length=nobj, dtype='f8')) # redshift
out.add_column(Column(name='CONTINUUM_COEFF', length=nobj, shape=(nssp_coeff,), dtype='f8'))
out.add_column(Column(name='CONTINUUM_RCHI2', length=nobj, dtype='f4')) # reduced chi2
Expand Down Expand Up @@ -2126,8 +2127,8 @@ def continuum_specfit(self, data):
bestfit = bestsspflux.dot(coeff)
meanage = self.get_meanage(coeff)

flam_ivar = np.hstack(data['ivar'])
dn4000, dn4000_ivar = self.get_dn4000(specwave, specflux, flam_ivar=flam_ivar, # specivar is line-masked!
flam_ivar = np.hstack(data['ivar']) # specivar is line-masked!
dn4000, dn4000_ivar = self.get_dn4000(specwave, specflux, flam_ivar=flam_ivar,
redshift=redshift, rest=False)
dn4000_model, _ = self.get_dn4000(specwave, bestfit, redshift=redshift, rest=False)

Expand Down Expand Up @@ -2284,6 +2285,7 @@ def continuum_joint_specfit(self, data):
# smooth continuum correction), which will guarantee that the aperture
# correction is always positive.
apcorr = data['phot']['nanomaggies'][1] / data['synthphot']['nanomaggies'][1]
data['apcorr'] = apcorr

# Prepare the reddened and unreddened SSP templates by redshifting and
# normalizing. Note that we ignore templates which are older than the
Expand Down Expand Up @@ -2442,24 +2444,28 @@ def continuum_joint_specfit(self, data):
self.log.info('Insufficient S/N to compute vdisp; adopting vdisp={:.2f} km/s'.format(self.vdisp_nominal))
vdispbest, vdispivar = self.vdisp_nominal, 0.0

pdb.set_trace()

# Get the final set of coefficients and chi2 at the best-fitting
# reddening and velocity dispersion.
bestsspflux, _ = self.SSP2data(self.sspflux[:, agekeep], self.sspwave, redshift=redshift,
specwave=data['wave'], specres=data['res'],
AV=AVbest, vdisp=vdispbest, cameras=data['cameras'],
south=data['photsys'] == 'S', synthphot=False)
# dev - synthesize photometry
bestsspflux, bestphot = self.SSP2data(self.sspflux[:, agekeep], self.sspwave, redshift=redshift,
specwave=data['wave'], specres=data['res'],
AV=AVbest, vdisp=vdispbest, cameras=data['cameras'],
south=data['photsys'] == 'S', synthphot=True)
bestsspflux = np.concatenate(bestsspflux, axis=0)
coeff, chi2min = self._call_nnls(bestsspflux, specflux, specivar)
chi2min /= np.sum(specivar > 0) # dof???
bestflam = bestphot['flam'].data*self.massnorm*self.fluxnorm

# dev - prepend the photometry
coeff, chi2min = self._call_nnls(np.vstack((bestflam, bestsspflux)),
np.hstack((objflam, specflux*apcorr)),
np.hstack((objflamivar, specivar/apcorr**2)))
chi2min /= (np.sum(objflamivar > 0) + np.sum(specivar > 0)) # dof???

# Get the light-weighted age and DN(4000).
bestfit = bestsspflux.dot(coeff)
meanage = self.get_meanage(coeff)

flam_ivar = np.hstack(data['ivar'])
dn4000, dn4000_ivar = self.get_dn4000(specwave, specflux, flam_ivar=flam_ivar, # specivar is line-masked!
flam_ivar = np.hstack(data['ivar']) # specivar is line-masked!
dn4000, dn4000_ivar = self.get_dn4000(specwave, specflux, flam_ivar=flam_ivar,
redshift=redshift, rest=False)
dn4000_model, _ = self.get_dn4000(specwave, bestfit, redshift=redshift, rest=False)

Expand Down Expand Up @@ -2511,7 +2517,7 @@ def continuum_joint_specfit(self, data):
if np.all(coeff == 0):
_smooth_continuum = np.zeros_like(bestfit)
else:
_smooth_continuum, _ = self.smooth_continuum(specwave, specflux - bestfit, specivar,
_smooth_continuum, _ = self.smooth_continuum(specwave, apcorr*specflux - bestfit, specivar/apcorr**2,
redshift, linemask=linemask, png=png)

# Unpack the continuum into individual cameras.
Expand All @@ -2522,24 +2528,25 @@ def continuum_joint_specfit(self, data):

## Like above, but with per-camera smoothing.
#smooth_continuum = self.smooth_residuals(
# continuummodel, data['wave'], data['flux'],
# continuummodel, data['wave'], apcorr*data['flux'],
# data['ivar'], data['linemask'], percamera=False)

if False:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(2, 1)
for icam in np.arange(len(data['cameras'])): # iterate over cameras
resid = data['flux'][icam]-continuummodel[icam]
resid = apcorr*data['flux'][icam]-continuummodel[icam]
ax[0].plot(data['wave'][icam], resid)
ax[1].plot(data['wave'][icam], resid-smooth_continuum[icam])
for icam in np.arange(len(data['cameras'])): # iterate over cameras
resid = data['flux'][icam]-continuummodel[icam]
resid = apcorr*data['flux'][icam]-continuummodel[icam]
pix_emlines = np.logical_not(data['linemask'][icam]) # affected by line = True
ax[0].scatter(data['wave'][icam][pix_emlines], resid[pix_emlines], s=30, color='red')
ax[0].plot(data['wave'][icam], smooth_continuum[icam], color='k', alpha=0.7, lw=2)
plt.savefig('junk.png')

# Pack it in and return.
result['APCORR'] = apcorr
result['CONTINUUM_COEFF'][0][0:nage] = coeff
result['CONTINUUM_RCHI2'][0] = chi2min
result['CONTINUUM_AV'][0] = AVbest # * u.mag
Expand Down
41 changes: 27 additions & 14 deletions py/fastspecfit/emlines.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -981,8 +981,8 @@ def fit(self, data, continuummodel, smooth_continuum, synthphot=True,
# best-fitting model (per-camera) below.
redshift = data['zredrock']
emlinewave = np.hstack(data['wave'])
oemlineivar = np.hstack(data['ivar'])
specflux = np.hstack(data['flux'])
oemlineivar = np.hstack(data['ivar'])/data['apcorr']**2
specflux = np.hstack(data['flux'])*data['apcorr']
resolution_matrix = data['res']
camerapix = data['camerapix']

Expand Down Expand Up @@ -1593,7 +1593,9 @@ def qa_fastspec(self, data, fastspec, metadata, coadd_type='healpix',
else:
pass

apcorr = fastspec['APCORR']
redshift = fastspec['CONTINUUM_Z']

stackwave = np.hstack(data['wave'])

# rebuild the best-fitting spectroscopic and photometric models
Expand All @@ -1604,14 +1606,23 @@ def qa_fastspec(self, data, fastspec, metadata, coadd_type='healpix',
vdisp=fastspec['CONTINUUM_VDISP'],
coeff=fastspec['CONTINUUM_COEFF'],
synthphot=False)

residuals = [data['flux'][icam] - continuum[icam] for icam in np.arange(len(data['cameras']))]
print('This bit is new!')
continuum_phot, synthmodelphot = self.SSP2data(
self.sspflux, self.sspwave, redshift=redshift,
synthphot=True, AV=fastspec['CONTINUUM_AV'], #test=True,
vdisp=fastspec['CONTINUUM_VDISP'],
coeff=fastspec['CONTINUUM_COEFF'] * self.massnorm)
continuum_wave_phot = self.sspwave * (1 + redshift)

residuals = [apcorr*data['flux'][icam] - continuum[icam] for icam in np.arange(len(data['cameras']))]

if np.all(fastspec['CONTINUUM_COEFF'] == 0):
_smooth_continuum = np.zeros_like(stackwave)
else:
_smooth_continuum, _ = self.smooth_continuum(np.hstack(data['wave']), np.hstack(residuals),
np.hstack(data['ivar']), redshift=redshift,
linemask=np.hstack(data['linemask']))

smooth_continuum = []
for campix in data['camerapix']:
smooth_continuum.append(_smooth_continuum[campix[0]:campix[1]])
Expand Down Expand Up @@ -1697,23 +1708,23 @@ def qa_fastspec(self, data, fastspec, metadata, coadd_type='healpix',
#ymin = np.zeros(len(data['cameras']))
#ymax = np.zeros(len(data['cameras']))
for ii in np.arange(len(data['cameras'])): # iterate over cameras
sigma, good = ivar2var(data['ivar'][ii], sigma=True, allmasked_ok=True)
sigma, good = ivar2var(data['ivar'][ii]/apcorr**2, sigma=True, allmasked_ok=True)

bigax1.fill_between(data['wave'][ii], data['flux'][ii]-sigma,
data['flux'][ii]+sigma, color=col1[ii])
bigax1.fill_between(data['wave'][ii], apcorr*data['flux'][ii]-sigma,
apcorr*data['flux'][ii]+sigma, color=col1[ii])
#bigax1.plot(data['wave'][ii], continuum_nodust[ii], alpha=0.5, color='k')
#bigax1.plot(data['wave'][ii], smooth_continuum[ii], color='gray')#col3[ii])#, alpha=0.3, lw=2)#, color='k')
#bigax1.plot(data['wave'][ii], continuum[ii], color=col2[ii])
bigax1.plot(data['wave'][ii], continuum[ii]+smooth_continuum[ii], color=col2[ii])

# get the robust range
filtflux = median_filter(data['flux'][ii], 51, mode='nearest')
filtflux = median_filter(apcorr*data['flux'][ii], 51, mode='nearest')
#filtflux = median_filter(data['flux'][ii] - _emlinemodel[ii], 51, mode='nearest')
#perc = np.percentile(filtflux[data['ivar'][ii] > 0], [5, 95])
#sigflux = np.std(data['flux'][ii][data['ivar'][ii] > 0])
#sigflux = np.std(apcorr*data['flux'][ii][data['ivar'][ii] > 0])
I = data['ivar'][ii] > 0
if np.sum(I) > 0:
sigflux = np.diff(np.percentile(data['flux'][ii][I], [25, 75]))[0] / 1.349 # robust
sigflux = np.diff(np.percentile(apcorr*data['flux'][ii][I], [25, 75]))[0] / 1.349 # robust
else:
sigflux = 0.0
#sigflux = np.std(filtflux[data['ivar'][ii] > 0])
Expand Down Expand Up @@ -1743,6 +1754,8 @@ def qa_fastspec(self, data, fastspec, metadata, coadd_type='healpix',
bigax1.plot(stackwave, _smooth_continuum, color='gray')#col3[ii])#, alpha=0.3, lw=2)#, color='k')
bigax1.plot(stackwave, np.hstack(continuum), color='k', alpha=0.1)#col3[ii])#, alpha=0.3, lw=2)#, color='k')

#bigax1.plot(continuum_wave_phot, continuum_phot / self.massnorm, color='orange', lw=4)

bigax1.text(0.03, 0.9, 'Observed Spectrum + Continuum Model',
ha='left', va='center', transform=bigax1.transAxes, fontsize=30)
if not self.nolegend:
Expand Down Expand Up @@ -1770,10 +1783,10 @@ def qa_fastspec(self, data, fastspec, metadata, coadd_type='healpix',
ymin, ymax = 1e6, -1e6
for ii in np.arange(len(data['cameras'])): # iterate over cameras
emlinewave = data['wave'][ii]
emlineflux = data['flux'][ii] - continuum[ii] - smooth_continuum[ii]
emlineflux = apcorr*data['flux'][ii] - continuum[ii] - smooth_continuum[ii]
emlinemodel = _emlinemodel[ii]

emlinesigma, good = ivar2var(data['ivar'][ii], sigma=True, allmasked_ok=True)
emlinesigma, good = ivar2var(data['ivar'][ii]/apcorr**2, sigma=True, allmasked_ok=True)
emlinewave = emlinewave[good]
emlineflux = emlineflux[good]
emlinesigma = emlinesigma[good]
Expand Down Expand Up @@ -1884,10 +1897,10 @@ def qa_fastspec(self, data, fastspec, metadata, coadd_type='healpix',
# iterate over cameras
for ii in np.arange(len(data['cameras'])): # iterate over cameras
emlinewave = data['wave'][ii]
emlineflux = data['flux'][ii] - continuum[ii] - smooth_continuum[ii]
emlineflux = apcorr*data['flux'][ii] - continuum[ii] - smooth_continuum[ii]
emlinemodel = _emlinemodel[ii]

emlinesigma, good = ivar2var(data['ivar'][ii], sigma=True, allmasked_ok=True)
emlinesigma, good = ivar2var(data['ivar'][ii]/apcorr**2, sigma=True, allmasked_ok=True)
emlinewave = emlinewave[good]
emlineflux = emlineflux[good]
emlinesigma = emlinesigma[good]
Expand Down
2 changes: 1 addition & 1 deletion py/fastspecfit/fastspecfit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -458,7 +458,7 @@ def build_webqa(CFit, EMFit, data, fastfit, metadata, coadd_type='healpix',
_emlinemodel_oneline.append(_emlinemodel_oneline1)

# Grab the viewer cutout.
width = int(60 / 0.262) # =1 arcmin
width = int(40 / 0.262) # =1 arcmin
height = int(width / 1.3) # 3:2 aspect ratio

cutoutpng = os.path.join('/tmp', 'tmp.'+os.path.basename(pngfile))
Expand Down

0 comments on commit 5f60414

Please sign in to comment.