Skip to content

Commit

Permalink
add --debug-plots option
Browse files Browse the repository at this point in the history
  • Loading branch information
@moustakas
moustakas committed Jul 31, 2023
1 parent f6a8dbc commit 58bad63
Show file tree
Hide file tree
Showing 3 changed files with 54 additions and 27 deletions.
50 changes: 31 additions & 19 deletions py/fastspecfit/continuum.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1430,7 +1430,7 @@ def templates2data(self, _templateflux, _templatewave, redshift=0.0, dluminosity

@staticmethod
def call_nnls(modelflux, flux, ivar, xparam=None, debug=False,
interpolate_coeff=False, xlabel=None, log=None):
interpolate_coeff=False, xlabel=None, png=None, log=None):
"""Wrapper on nnls.
Works with both spectroscopic and photometric input and with both 2D and
Expand Down Expand Up @@ -1480,11 +1480,14 @@ def call_nnls(modelflux, flux, ivar, xparam=None, debug=False,

try:
imin = find_minima(chi2grid)[0]
xbest, xerr, chi2min, warn = minfit(xparam[imin-1:imin+2], chi2grid[imin-1:imin+2])
if debug:
xbest, xerr, chi2min, warn, (a, b, c) = minfit(xparam[imin-1:imin+2], chi2grid[imin-1:imin+2], return_coeff=True)
else:
xbest, xerr, chi2min, warn = minfit(xparam[imin-1:imin+2], chi2grid[imin-1:imin+2])
except:
errmsg = 'A problem was encountered minimizing chi2.'
log.warning(errmsg)
imin, xbest, xerr, chi2min, warn = 0, 0.0, 0.0, 0.0, 1
imin, xbest, xerr, chi2min, warn, (a, b, c) = 0, 0.0, 0.0, 0.0, 1, (0., 0., 0.)

if warn == 0:
xivar = 1.0 / xerr**2
Expand All @@ -1510,29 +1513,37 @@ def call_nnls(modelflux, flux, ivar, xparam=None, debug=False,

if debug:
if xivar > 0:
leg = r'${:.3f}\pm{:.3f}\ (\chi^2_{{min}}={:.2f})$'.format(xbest, 1/np.sqrt(xivar), chi2min)
leg = r'${:.1f}\pm{:.1f}$'.format(xbest, 1 / np.sqrt(xivar))
#leg = r'${:.3f}\pm{:.3f}\ (\chi^2_{{min}}={:.2f})$'.format(xbest, 1/np.sqrt(xivar), chi2min)
else:
leg = r'${:.3f}$'.format(xbest)

if xlabel:
leg = f'{xlabel} = {leg}'

import matplotlib.pyplot as plt
import seaborn as sns
sns.set(context='talk', style='ticks', font_scale=0.8)
fig, ax = plt.subplots()
ax.scatter(xparam, chi2grid)
ax.scatter(xparam[imin-1:imin+2], chi2grid[imin-1:imin+2], color='red')
#ax.set_ylim(chi2min*0.95, np.max(chi2grid[imin-1:imin+2])*1.05)
#ax.plot(xx, np.polyval([aa, bb, cc], xx), ls='--')
ax.axvline(x=xbest, color='k')
if xivar > 0:
ax.axhline(y=chi2min, color='k')
#ax.set_yscale('log')
#ax.set_ylim(chi2min, 63.3)
ax.scatter(xparam, chi2grid-chi2min, marker='s', s=50, color='gray', edgecolor='k')
#ax.scatter(xparam[imin-1:imin+2], chi2grid[imin-1:imin+2]-chi2min, color='red')
if not np.all(np.array([a, b, c]) == 0):
ax.plot(xparam, np.polyval([a, b, c], xparam)-chi2min, lw=2, ls='--')
#ax.axvline(x=xbest, color='k')
#if xivar > 0:
# ax.axhline(y=chi2min, color='k')
if xlabel:
ax.set_xlabel(xlabel)
#ax.text(0.03, 0.9, '{}={}'.format(xlabel, leg), ha='left',
# va='center', transform=ax.transAxes)
ax.text(0.03, 0.9, leg, ha='left', va='center', transform=ax.transAxes)
ax.set_ylabel(r'$\chi^2$')
ax.text(0.9, 0.9, leg, ha='right', va='center', transform=ax.transAxes)
ax.set_ylabel(r'$\Delta\chi^2$')
#ax.set_ylabel(r'$\chi^2 - {:.1f}$'.format(chi2min))
#fig.savefig('qa-chi2min.png')
fig.savefig('desi-users/ioannis/tmp/qa-chi2min.png')
fig.tight_layout()
if png:
log.info(f'Writing {png}')
fig.savefig(png)

return chi2min, xbest, xivar, bestcoeff

Expand Down Expand Up @@ -1742,7 +1753,7 @@ def _kcorr_and_absmag(filters_out, band_shift):
return kcorr, absmag, ivarabsmag, bestmaggies, lums, cfluxes

def continuum_specfit(data, result, templatecache, nophoto=False, constrain_age=False,
fastphot=False, log=None, verbose=False):
fastphot=False, log=None, verbose=False, debug_plots=False):
"""Fit the non-negative stellar continuum of a single spectrum.
Parameters
Expand Down Expand Up @@ -1905,9 +1916,10 @@ def _younger_than_universe(age, tuniv, agepad=0.5):
vdispchi2min, vdispbest, vdispivar, _ = CTools.call_nnls(
ztemplateflux_vdisp[Ivdisp, :, :],
specflux[Ivdisp], specivar[Ivdisp],
xparam=templatecache['vdisp'], xlabel=r'$\sigma$ (km/s)', debug=False, log=log)
xparam=templatecache['vdisp'], xlabel=r'$\sigma$ (km/s)', log=log,
debug=debug_plots, png='deltachi2-vdisp.png')
log.info('Fitting for the velocity dispersion took {:.2f} seconds.'.format(time.time()-t0))

if vdispivar > 0:
# Require vdisp to be measured with S/N>1, which protects
# against tiny ivar becomming infinite in the output table.
Expand Down
8 changes: 5 additions & 3 deletions py/fastspecfit/fastspecfit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,8 @@ def _assign_units_to_columns(fastfit, metadata, Spec, templates, fastphot, stack

def fastspec_one(iobj, data, out, meta, templates, log=None,
minspecwave=3500., maxspecwave=9900., broadlinefit=True,
fastphot=False, stackfit=False, nophoto=False, percamera_models=False):
fastphot=False, stackfit=False, nophoto=False,
percamera_models=False, debug_plots=False):
"""Multiprocessing wrapper to run :func:`fastspec` on a single object.
"""
Expand All @@ -60,7 +61,7 @@ def fastspec_one(iobj, data, out, meta, templates, log=None,

continuummodel, smooth_continuum = continuum_specfit(data, out, templatecache,
fastphot=fastphot, nophoto=nophoto,
log=log)
debug_plots=debug_plots, log=log)

# Optionally fit the emission-line spectrum.
if fastphot:
Expand Down Expand Up @@ -125,6 +126,7 @@ def parse(options=None, log=None):
parser.add_argument('--mapdir', type=str, default=None, help='Optional directory name for the dust maps.')
parser.add_argument('--dr9dir', type=str, default=None, help='Optional directory name for the DR9 photometry.')
parser.add_argument('--specproddir', type=str, default=None, help='Optional directory name for the spectroscopic production.')
parser.add_argument('--debug-plots', action='store_true', help='Generate a variety of debugging plots (written to $PWD).')
parser.add_argument('--verbose', action='store_true', help='Be verbose (for debugging purposes).')

if log is None:
Expand Down Expand Up @@ -225,7 +227,7 @@ def fastspec(fastphot=False, stackfit=False, args=None, comm=None, verbose=False
t0 = time.time()
fitargs = [(iobj, data[iobj], out[iobj], meta[iobj], templates, log,
minspecwave, maxspecwave, args.broadlinefit, fastphot, stackfit,
args.nophoto, args.percamera_models)
args.nophoto, args.percamera_models, args.debug_plots)
for iobj in np.arange(Spec.ntargets)]
if args.mp > 1:
import multiprocessing
Expand Down
23 changes: 18 additions & 5 deletions py/fastspecfit/util.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -231,7 +231,7 @@ def find_minima(x):
return ii[jj]


def minfit(x, y):
def minfit(x, y, return_coeff=False):
"""Fits y = y0 + ((x-x0)/xerr)**2
See redrock.zwarning.ZWarningMask.BAD_MINFIT for zwarn failure flags
Expand All @@ -244,17 +244,27 @@ def minfit(x, y):
(tuple): (x0, xerr, y0, zwarn) where zwarn=0 is good fit.
"""
a, b, c = 0., 0., 0.
if len(x) < 3:
return (-1,-1,-1,ZWarningMask.BAD_MINFIT)
if return_coeff:
return (-1,-1,-1,ZWarningMask.BAD_MINFIT,(a,b,c))
else:
return (-1,-1,-1,ZWarningMask.BAD_MINFIT)

try:
#- y = a x^2 + b x + c
a,b,c = np.polyfit(x,y,2)
except np.linalg.LinAlgError:
return (-1,-1,-1,ZWarningMask.BAD_MINFIT)
if return_coeff:
return (-1,-1,-1,ZWarningMask.BAD_MINFIT,(a,b,c))
else:
return (-1,-1,-1,ZWarningMask.BAD_MINFIT)

if a == 0.0:
return (-1,-1,-1,ZWarningMask.BAD_MINFIT)
if return_coeff:
return (-1,-1,-1,ZWarningMask.BAD_MINFIT,(a,b,c))
else:
return (-1,-1,-1,ZWarningMask.BAD_MINFIT)

#- recast as y = y0 + ((x-x0)/xerr)^2
x0 = -b / (2*a)
Expand All @@ -272,7 +282,10 @@ def minfit(x, y):
xerr = 1 / np.sqrt(-a)
zwarn |= ZWarningMask.BAD_MINFIT

return (x0, xerr, y0, zwarn)
if return_coeff:
return (x0, xerr, y0, zwarn,(a,b,c))
else:
return (x0, xerr, y0, zwarn)

class TabulatedDESI(object):
"""
Expand Down

0 comments on commit 58bad63

Please sign in to comment.