From f26690421452a5d6a9234a8b5aa19133460a1235 Mon Sep 17 00:00:00 2001 From: jdbuhler Date: Fri, 30 Aug 2024 23:19:48 -0400 Subject: [PATCH 01/20] Reimplement the DESI resolution matrix using Numba for the multiply, which yields a nice speedup. Merge the needed functionality of the emline_fit sparse matrix, so that we don't have to pass two sepearate resolution matrices to different parts of the code. --- py/fastspecfit/continuum.py | 8 +- py/fastspecfit/emline_fit/__init__.py | 2 - py/fastspecfit/emline_fit/interface.py | 10 +- py/fastspecfit/emline_fit/sparse_rep.py | 89 +--------------- py/fastspecfit/emlines.py | 2 +- py/fastspecfit/io.py | 17 ++- py/fastspecfit/resolution.py | 133 ++++++++++++++++++++++++ 7 files changed, 152 insertions(+), 109 deletions(-) create mode 100644 py/fastspecfit/resolution.py diff --git a/py/fastspecfit/continuum.py b/py/fastspecfit/continuum.py index 7ab58838..95667f38 100644 --- a/py/fastspecfit/continuum.py +++ b/py/fastspecfit/continuum.py @@ -13,8 +13,8 @@ from fastspecfit.photometry import Photometry from fastspecfit.templates import Templates from fastspecfit.util import (C_LIGHT, FLUXNORM, - trapz_rebin, trapz_rebin_pre, quantile, - median, sigmaclip) + trapz_rebin, trapz_rebin_pre, + quantile, median, sigmaclip) class ContinuumTools(object): @@ -756,13 +756,13 @@ def continuum_to_spectroscopy(self, contmodel): specres = self.data['res'] modelflux = np.empty(self.wavelen) - + for icam, (s, e) in enumerate(camerapix): resampflux = trapz_rebin(self.ztemplatewave, contmodel, specwave[icam], pre=self.spec_pre[icam]) - modelflux[s:e] = specres[icam] @ resampflux + specres[icam].dot(resampflux, out=modelflux[s:e]) return modelflux diff --git a/py/fastspecfit/emline_fit/__init__.py b/py/fastspecfit/emline_fit/__init__.py index 9ac3417b..448de7aa 100644 --- a/py/fastspecfit/emline_fit/__init__.py +++ b/py/fastspecfit/emline_fit/__init__.py @@ -6,5 +6,3 @@ ) from .params_mapping import ParamsMapping as EMLine_ParamsMapping - -from .sparse_rep import ResMatrix as EMLine_Resolution diff --git a/py/fastspecfit/emline_fit/interface.py b/py/fastspecfit/emline_fit/interface.py index fda86691..1a334dd4 100644 --- a/py/fastspecfit/emline_fit/interface.py +++ b/py/fastspecfit/emline_fit/interface.py @@ -164,14 +164,14 @@ def jacobian(self, free_parameters): ibw, self.redshift, self.line_wavelengths, - self.resolution_matrices[icam].ndiag()) + self.resolution_matrices[icam].ndiag) # ignore any columns corresponding to fixed parameters endpts = idealJac[0] endpts[self.params_mapping.fixedMask(), :] = 0 jacs.append( mulWMJ(self.obs_weights[s:e], - self.resolution_matrices[icam].data, + self.resolution_matrices[icam].rowdata(), idealJac) ) nBins = np.sum(np.diff(self.camerapix)) @@ -445,7 +445,7 @@ def _build_model_core(line_parameters, # convolve model with resolution matrix and store in # this camera's subrange of model_fluxes - resolution_matrices[icam].matvec(mf, model_fluxes[s:e]) + resolution_matrices[icam].dot(mf, model_fluxes[s:e]) # @@ -494,11 +494,11 @@ def _build_multimodel_core(line_parameters, log_obs_bin_edges[s+icam:e+icam+1], redshift, ibw, - resolution_matrices[icam].ndiag()) + resolution_matrices[icam].ndiag) # convolve each line's waveform with resolution matrix endpts, M = mulWMJ(np.ones(e - s), - resolution_matrices[icam].data, + resolution_matrices[icam].rowdata(), line_models) # adjust endpoints to reflect camera range diff --git a/py/fastspecfit/emline_fit/sparse_rep.py b/py/fastspecfit/emline_fit/sparse_rep.py index 2936adef..bcdea032 100644 --- a/py/fastspecfit/emline_fit/sparse_rep.py +++ b/py/fastspecfit/emline_fit/sparse_rep.py @@ -1,7 +1,6 @@ # -# Sparse representations of resolution matrices and -# Jacobian of objective for emission line fitting -# (Ideal Jacobian rep is generated in EMLines_jacobian()) +# Sparse representations of Jacobian of objective for emission line +# fitting (Ideal Jacobian rep is generated in EMLines_jacobian()) # import numpy as np @@ -11,90 +10,6 @@ from .params_mapping import ParamsMapping - -# -# resolution matrix -# A resolution matrix M of size nrow x nrow is stored as a 2D array A of -# size nrow x ndiag, where ndiag is the number of nonzero diagonals -# (which must be odd). The rows of A are M's rows, but with only the -# nonzero entries stored. The nonzero entries on row i run from -# j = i - diag//2 to i + diag//2, so -# M[i,j] = A[i, j - (i - diag//2)] -# -class ResMatrix(object): - - def __init__(self, D): - self.data = self._from_dia_matrix(D) - - def ndiag(self): - return self.data.shape[1] - - def matvec(self, v, w): - self._matvec(self.data, v, w) - - # for compatibility - def dot(self, v): - w = np.empty(self.data.shape[0]) - self.matvec(v, w) - return w - - # - # _from_dia_matrix() - # Convert a diagonally sparse matrix M in the form - # stored by DESI into a sparse row rerpesentation. - # - # Input M is represented as a 2D array D of size ndiag x nrow, - # whose rows are M's diagonals: - # M[i,j] = D[ndiag//2 - (j - i), j] - # ndiag is assumed to be odd, and entries in D that would be - # outside the bounds of M are ignored. - # - @staticmethod - @jit(nopython=True, fastmath=False, nogil=True) - def _from_dia_matrix(D): - ndiag, nrow = D.shape - hdiag = ndiag//2 - - A = np.empty((nrow, ndiag), dtype=D.dtype) - - for i in range(nrow): - # min and max column for row - jmin = np.maximum(i - hdiag, 0) - jmax = np.minimum(i + hdiag, nrow - 1) - - for j in range(jmin, jmax + 1): - A[i, j - i + hdiag] = D[hdiag + i - j, j] - - return A - - # - # _matvec() - # Compute the matrix-vector product M * v, where - # M is a row-sparse matrix with a limited - # number of diagonals created by - # dia_to_row_matrix(). - # - # w is an output parameter - # - @staticmethod - @jit(nopython=True, fastmath=False, nogil=True) - def _matvec(M, v, w): - nrow, ndiag = M.shape - hdiag = ndiag//2 - - for i in range(nrow): - jmin = np.maximum(i - hdiag, 0) - jmax = np.minimum(i + hdiag, nrow - 1) - - acc = 0. - for j in range(jmin, jmax + 1): - acc += M[i, j - i + hdiag] * v[j] - - w[i] = acc - - -################################################################# - # # Sparse Jacobian of objective function. For # each camera's pixel range, the Jacobian diff --git a/py/fastspecfit/emlines.py b/py/fastspecfit/emlines.py index 96658cf9..d0794982 100644 --- a/py/fastspecfit/emlines.py +++ b/py/fastspecfit/emlines.py @@ -1100,7 +1100,7 @@ def emline_specfit(data, result, continuummodel, smooth_continuum, redshift = data['redshift'] camerapix = data['camerapix'] - resolution_matrix = data['res_emline'] + resolution_matrix = data['res'] # Combine pixels across all cameras emlinewave = np.hstack(data['wave']) diff --git a/py/fastspecfit/io.py b/py/fastspecfit/io.py index ddcf4d80..b32c7bab 100644 --- a/py/fastspecfit/io.py +++ b/py/fastspecfit/io.py @@ -637,7 +637,7 @@ def read(self, mp_pool, fastphot=False, synthphot=True, Three-element list of `numpy.ndarray` inverse variance spectra, one for each camera. res : :class:`list` - Three-element list of :class:`desispec.resolution.Resolution` + Three-element list of :class:`fastspecfit.resolution.Resolution` objects, one for each camera. snr : `numpy.ndarray` Median per-pixel signal-to-noise ratio in the grz cameras. @@ -681,8 +681,7 @@ def read(self, mp_pool, fastphot=False, synthphot=True, from desispec.coaddition import coadd_cameras from desispec.io import read_spectra from desiutil.dust import SFDMap - from desispec.resolution import Resolution - from fastspecfit.emline_fit import EMLine_Resolution + from fastspecfit.resolution import Resolution t0 = time.time() @@ -782,7 +781,7 @@ def read(self, mp_pool, fastphot=False, synthphot=True, 'coadd_flux': coadd_spec.flux[coadd_cams][iobj, :], 'coadd_ivar': coadd_spec.ivar[coadd_cams][iobj, :], 'coadd_res': [Resolution(coadd_spec.resolution_data[coadd_cams][iobj, :])], - 'coadd_res_emline': [EMLine_Resolution(coadd_spec.resolution_data[coadd_cams][iobj, :])], + 'coadd_res_emline': [Resolution(coadd_spec.resolution_data[coadd_cams][iobj, :])], } mpargs.append({ @@ -821,8 +820,7 @@ def one_spectrum(iobj, specdata, meta, ebv, fastphot, extinction. Also flag pixels which may be affected by emission lines. """ - from desispec.resolution import Resolution - from fastspecfit.emline_fit import EMLine_Resolution + from fastspecfit.resolution import Resolution from fastspecfit.util import mwdust_transmission, median from fastspecfit.linemasker import LineMasker @@ -949,8 +947,7 @@ def one_spectrum(iobj, specdata, meta, ebv, fastphot, specdata['mask'].append(specdata['mask0'][icam]) specdata['res'].append(Resolution(res)) - specdata['res_emline'].append(EMLine_Resolution(res)) - + if len(cameras) == 0: errmsg = 'No good data, which should never happen.' log.critical(errmsg) @@ -1058,7 +1055,7 @@ def read_stacked(self, stackfiles, firsttarget=0, ntargets=None, Three-element list of `numpy.ndarray` inverse variance spectra, one for each camera. res : :class:`list` - Three-element list of :class:`desispec.resolution.Resolution` + Three-element list of :class:`fastspecfit.resolution.Resolution` objects, one for each camera. snr : `numpy.ndarray` Median per-pixel signal-to-noise ratio in the grz cameras. @@ -1098,7 +1095,7 @@ def read_stacked(self, stackfiles, firsttarget=0, ntargets=None, """ from astropy.table import vstack - from desispec.resolution import Resolution + from fastspecfit.resolution import Resolution if stackfiles is None: errmsg = 'At least one stackfiles file is required.' diff --git a/py/fastspecfit/resolution.py b/py/fastspecfit/resolution.py new file mode 100644 index 00000000..a140986d --- /dev/null +++ b/py/fastspecfit/resolution.py @@ -0,0 +1,133 @@ +import numpy as np +from numba import jit + +class Resolution(object): + """Resolution matrix, in the style of desispec.resolution.Resolution + + The base implementation is analogous to scipy.sparse.dia_matrix, + storing a [ndiag x dim] 2D array giving the values on diagonals + [ndiag//2 ... -ndiag//2] of a dim x dim matrix. We reimplement + the dot() method in Numba, special-casing for this specific + pattern of diagonal offsets; the result is faster than using + dia_matrix, even though the latter's dot method is in C. + + We also provide conversion to a 2D array giving the nonzero + entries in each *row* of the matrix, which is used by emline + fitting code. + + """ + + def __init__(self, data): + """ + Parameters + ---------- + data : :class:`np.ndarray` [ndiag x dim] + Array of values for diagonals [ndiag//2 .. -ndiag//2]; + note that some entries at the ends of rows are ignored. + All other diagonals are assumed to be zero. + + """ + ndiag, dim = data.shape + + self.shape = (dim, dim) + self.ndiag = ndiag + self.data = data + + self.rows = None # compute on first use + + + def rowdata(self): + """ + Compute sparse row matrix from diagonal representation. + + Returns + ------- + 2D array of size [dim x ndiag] with nonzero entries + for each row of the matrix. + + """ + if self.rows is None: + self.rows = self._dia_to_rows(self.data) + return self.rows + + + def dot(self, v, out=None): + """ + Compute our dot product with vector v, storing + result in out (which is created if None). + + Parameters + ---------- + v : :class:`np.ndarray` [dim] + vector of values to dot with us. + out : :class:`np.ndarray` [dim] or None + array to hold output; created if None. + + Returns + ------- + array of length [dim] holding output. + + """ + if out is None: + out = np.empty(self.shape[0], dtype=v.dtype) + + return self._matvec(self.data, v, out) + + + @staticmethod + @jit(nopython=True, fastmath=False, nogil=True) + def _matvec(D, v, out): + """ + Compute matrix-vector product of a Resolution matrix A and a vector v. + + A has shape dim x dim, and its nonzero entries are given by + a matrix D of size ndiag x dim, providing values on diagonals + [ndiag//2 .. -ndiag//2], inclusive. v is an array of length dim. + + Return result in array out of size dim. + + """ + out[:] = 0. + + ndiag, dim = D.shape + for i in range(ndiag): + k = ndiag//2 - i # diagonal offset + + i_start = np.maximum(0, -k) + j_start = np.maximum(0, k) + j_end = np.minimum(dim + k, dim) + + for n in range(j_end - j_start): + out[i_start + n] += D[i, j_start + n] * v[j_start + n] + + return out + + + @staticmethod + @jit(nopython=True, fastmath=False, nogil=True) + def _dia_to_rows(D): + """ + Convert a diagonally sparse matrix M in the form + stored by DESI into a sparse row representation. + + Input M is represented as a 2D array D of size ndiag x nrow, + whose rows are M's diagonals: + M[i,j] = D[ndiag//2 - (j - i), j] + ndiag is assumed to be odd, and entries in D that would be + outside the bounds of M are ignored. + """ + + ndiag, dim = D.shape + hdiag = ndiag//2 + + A = np.empty((dim, ndiag), dtype=D.dtype) + + for i in range(dim): + # min and max column for row + jmin = np.maximum(i - hdiag, 0) + jmax = np.minimum(i + hdiag, dim - 1) + + for j in range(jmin, jmax + 1): + A[i, j - (i - hdiag)] = D[hdiag + i - j, j] + + return A From 07b4c41daf11499abc93eb65732ebdbb5cce2c80 Mon Sep 17 00:00:00 2001 From: Jeremy Buhler Date: Sun, 1 Sep 2024 18:49:25 -0500 Subject: [PATCH 02/20] * clean up some stray garbage from emlines version of resolution matrices * don't allocate unused local reference to pixkms_bounds in ContinuumTools * emit a log message when we use an acceelerated FFT library * store frequently used, frozen lists in input data dictionary as fixed-sized tuples or arrays for efficiency --- py/fastspecfit/continuum.py | 7 +++--- py/fastspecfit/emlines.py | 4 ++-- py/fastspecfit/io.py | 46 +++++++++++++++++++----------------- py/fastspecfit/linemasker.py | 2 +- py/fastspecfit/templates.py | 3 ++- 5 files changed, 32 insertions(+), 30 deletions(-) diff --git a/py/fastspecfit/continuum.py b/py/fastspecfit/continuum.py index 95667f38..12b2f63f 100644 --- a/py/fastspecfit/continuum.py +++ b/py/fastspecfit/continuum.py @@ -32,7 +32,6 @@ def __init__(self, data, templates, phot, igm, ebv_guess=0.05, self.massnorm = 1e10 # stellar mass normalization factor [Msun] self.ebv_guess = 0.05 # [mag] - self.pixkms_bounds = templates.pixkms_bounds self.lg_atten = np.log(10.) * (-0.4 * templates.dust_klambda) # Cache the redshift-dependent factors (incl. IGM attenuation), @@ -341,7 +340,7 @@ def get_cflux(cwave): def templates2data(self, templateflux, templatewave, redshift=0., dluminosity=None, - vdisp=None, cameras=['b','r','z'], specwave=None, specres=None, + vdisp=None, cameras=np.array(['b','r','z']), specwave=None, specres=None, specmask=None, coeff=None, photsys=None, synthphot=True, stack_cameras=False, flamphot=False, debug=False, get_abmag=False): """Deprecated. Work-horse method for turning input templates into data-like @@ -748,13 +747,13 @@ def continuum_to_spectroscopy(self, contmodel): ------- modelflux : :class:`numpy.ndarray` [nwave] Observed-frame model spectrum at the instrumental resolution and - wavelength sampling given by `specwave` and `specres`. + wavelength sampling given by `specres` and `specwave`. """ camerapix = self.data['camerapix'] specwave = self.data['wave'] specres = self.data['res'] - + modelflux = np.empty(self.wavelen) for icam, (s, e) in enumerate(camerapix): diff --git a/py/fastspecfit/emlines.py b/py/fastspecfit/emlines.py index d0794982..3ca13057 100644 --- a/py/fastspecfit/emlines.py +++ b/py/fastspecfit/emlines.py @@ -576,8 +576,8 @@ def optimize(self, linemodel, obs_weights, redshift, line_wavelengths, - tuple(resolution_matrices), - camerapix, # for more efficient iteration + resolution_matrices, + camerapix, params_mapping, continuum_patches=continuum_patches) diff --git a/py/fastspecfit/io.py b/py/fastspecfit/io.py index b32c7bab..904a284e 100644 --- a/py/fastspecfit/io.py +++ b/py/fastspecfit/io.py @@ -781,8 +781,7 @@ def read(self, mp_pool, fastphot=False, synthphot=True, 'coadd_flux': coadd_spec.flux[coadd_cams][iobj, :], 'coadd_ivar': coadd_spec.ivar[coadd_cams][iobj, :], 'coadd_res': [Resolution(coadd_spec.resolution_data[coadd_cams][iobj, :])], - 'coadd_res_emline': [Resolution(coadd_spec.resolution_data[coadd_cams][iobj, :])], - } + } mpargs.append({ 'iobj': iobj, @@ -894,7 +893,6 @@ def one_spectrum(iobj, specdata, meta, ebv, fastphot, 'ivar': [], 'mask': [], 'res': [], - 'res_emline': [], 'snr': np.zeros(len(np.atleast_1d(specdata['cameras'])), 'f4'), 'linemask': [], 'linepix': [], @@ -953,27 +951,29 @@ def one_spectrum(iobj, specdata, meta, ebv, fastphot, log.critical(errmsg) raise ValueError(errmsg) - # clean up the data dictionary - for key in ('wave0', 'flux0', 'ivar0', 'mask0', 'res0'): - del specdata[key] - + # clean up unused items in data dictionary and + # freeze lists that will not be further modified + for key in ('wave', 'flux', 'ivar', 'mask', 'res'): + del specdata[key + '0'] + specdata[key] = tuple(specdata[key]) + # Pre-compute some convenience variables for "un-hstacking" # an "hstacked" spectrum. - specdata['cameras'] = cameras - specdata['npixpercamera'] = npixpercamera + specdata['cameras'] = np.array(cameras) + specdata['npixpercamera'] = np.array(npixpercamera) ncam = len(specdata['cameras']) - c_ends = np.cumsum(npixpercamera) - c_starts = c_ends - npixpercamera + c_ends = np.cumsum(specdata['npixpercamera']) + c_starts = c_ends - specdata['npixpercamera'] specdata['camerapix'] = np.zeros((ncam, 2), np.int16) specdata['camerapix'][:, 0] = c_starts specdata['camerapix'][:, 1] = c_ends - + # use the coadded spectrum to build a robust emission-line mask LM = LineMasker(emline_table) pix = LM.build_linemask( specdata['coadd_wave'], specdata['coadd_flux'], - specdata['coadd_ivar'], specdata['coadd_res_emline'], + specdata['coadd_ivar'], specdata['coadd_res'], uniqueid=specdata['uniqueid'], redshift=specdata['redshift'], debug_plots=debug_plots) @@ -982,7 +982,7 @@ def one_spectrum(iobj, specdata, meta, ebv, fastphot, # way to do it? for icam in range(ncam): camlinepix = {} - camlinemask = np.zeros(npixpercamera[icam], bool) + camlinemask = np.zeros(specdata['npixpercamera'][icam], bool) for linename in pix['coadd_linepix']: linepix = pix['coadd_linepix'][linename] # if the line is entirely off this camera, skip it @@ -1346,22 +1346,24 @@ def one_stacked_spectrum(iobj, specdata, synthphot): log.critical(errmsg) raise ValueError(errmsg) + # clean up unused items in data dictionary and + # freeze lists that will not be further modified + for key in ('wave', 'flux', 'ivar', 'mask', 'res'): + del specdata[key + '0'] + specdata[key] = tuple(specdata[key]) + # Pre-compute some convenience variables for "un-hstacking" # an "hstacked" spectrum. - specdata['cameras'] = cameras - specdata['npixpercamera'] = npixpercamera + specdata['cameras'] = np.array(cameras) + specdata['npixpercamera'] = np.array(npixpercamera) ncam = len(specdata['cameras']) - c_ends = np.cumsum(npixpercamera) - c_starts = c_ends - npixpercamera + c_ends = np.cumsum(specdata['npixpercamera']) + c_starts = c_ends - specdata['npixpercamera'] specdata['camerapix'] = np.zeros((ncam, 2), np.int16) specdata['camerapix'][:, 0] = c_starts specdata['camerapix'][:, 1] = c_ends - # clean up the data dictionary - for key in ('wave0', 'flux0', 'ivar0', 'mask0', 'res0'): - del specdata[key] - LM = LineMasker(phot, emline_table) coadd_linemask_dict = LM.build_linemask(specdata['coadd_wave'], specdata['coadd_flux'], specdata['coadd_ivar'], redshift=specdata['redshift']) diff --git a/py/fastspecfit/linemasker.py b/py/fastspecfit/linemasker.py index 53a542d4..14bea576 100644 --- a/py/fastspecfit/linemasker.py +++ b/py/fastspecfit/linemasker.py @@ -209,7 +209,7 @@ def build_linemask(self, wave, flux, ivar, resolution_matrix, redshift=0., ivar : :class:`numpy.ndarray` [npix] Inverse variance spectrum corresponding to `flux`. resolution_matrix : :class:`list` - List of :class:`fastspecfit.emline_fit.sparse_rep.ResMatrix` + List of :class:`fastspecfit.resolution.Resolution` objects, one per camera. redshift : :class:`float` Object redshift. diff --git a/py/fastspecfit/templates.py b/py/fastspecfit/templates.py index 4c708f6e..cf461de8 100644 --- a/py/fastspecfit/templates.py +++ b/py/fastspecfit/templates.py @@ -206,8 +206,9 @@ def init_ffts(self): if find_spec("mkl_fft") is not None: import mkl_fft._scipy_fft_backend as be sc_fft.set_global_backend(be) - + self.convolve = sc_sig.convolve + log.info('Using mkl_fft library for FFTs') else: self.convolve = sc_sig.oaconvolve From 7f36418d4f3be8f51e19147e54160301408bda3a Mon Sep 17 00:00:00 2001 From: Jeremy Buhler Date: Sun, 1 Sep 2024 19:55:22 -0500 Subject: [PATCH 03/20] * tiny cleanup of spectroscopy preprocessing --- py/fastspecfit/continuum.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/py/fastspecfit/continuum.py b/py/fastspecfit/continuum.py index 12b2f63f..0ffbbe5f 100644 --- a/py/fastspecfit/continuum.py +++ b/py/fastspecfit/continuum.py @@ -70,9 +70,7 @@ def __init__(self, data, templates, phot, igm, ebv_guess=0.05, self.wavelen = np.sum([len(w) for w in self.data['wave']]) # get preprocessing data to accelerate resampling - rspre = [] - for wave in self.data['wave']: - rspre.append(trapz_rebin_pre(wave)) + rspre = [ trapz_rebin_pre(w) for w in self.data['wave'] ] self.spec_pre = tuple(rspre) From 7ced18f6ece1b3684938f39ee429d5341d55dc94 Mon Sep 17 00:00:00 2001 From: Jeremy Buhler Date: Sun, 1 Sep 2024 20:47:34 -0500 Subject: [PATCH 04/20] minor cleanup to skip needless array copy in continuum objective --- py/fastspecfit/continuum.py | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/py/fastspecfit/continuum.py b/py/fastspecfit/continuum.py index 0ffbbe5f..b7b9abc0 100644 --- a/py/fastspecfit/continuum.py +++ b/py/fastspecfit/continuum.py @@ -831,8 +831,6 @@ def _stellar_objective(self, params, templateflux, vdisp = None templatecoeff = params[1:] - # FIXME: write results directly into resid array instead of - # allocating new fullmodel = self.build_stellar_continuum( templateflux, templatecoeff, ebv=ebv, vdisp=vdisp, conv_pre=conv_pre, @@ -857,17 +855,15 @@ def _stellar_objective(self, params, templateflux, if synthspec: modelflux = self.continuum_to_spectroscopy(fullmodel) spec_resid = resid[:resid_split] - spec_resid[:] = specflux - spec_resid -= modelflux - spec_resid *= specistd - + np.subtract(modelflux, specflux, spec_resid) + spec_resid *= specistd + if synthphot: modelflam = self.continuum_to_photometry(fullmodel) phot_resid = resid[resid_split:] - phot_resid[:] = objflam - phot_resid -= modelflam - phot_resid *= objflamistd - + np.subtract(modelflam, objflam, phot_resid) + phot_resid *= objflamistd + return resid From ebf455618d0fd5171ce2e84d64188dc9b1e486d0 Mon Sep 17 00:00:00 2001 From: Jeremy Buhler Date: Sun, 1 Sep 2024 22:21:54 -0500 Subject: [PATCH 05/20] remove unneeded import --- py/fastspecfit/util.py | 1 - 1 file changed, 1 deletion(-) diff --git a/py/fastspecfit/util.py b/py/fastspecfit/util.py index a93fa3a3..e273d3c2 100644 --- a/py/fastspecfit/util.py +++ b/py/fastspecfit/util.py @@ -5,7 +5,6 @@ General utilities. """ -import os import numpy as np import numba From 14ef0b12924697db20b04d532feb7bcced237315 Mon Sep 17 00:00:00 2001 From: jdbuhler Date: Tue, 3 Sep 2024 10:36:53 -0400 Subject: [PATCH 06/20] * fix missed uses of old res_emline object in QA * update documentation of sparse emline Jacobian to std --- py/fastspecfit/emline_fit/sparse_rep.py | 233 +++++++++++++----------- py/fastspecfit/qa.py | 4 +- 2 files changed, 131 insertions(+), 106 deletions(-) diff --git a/py/fastspecfit/emline_fit/sparse_rep.py b/py/fastspecfit/emline_fit/sparse_rep.py index bcdea032..02a8ec3b 100644 --- a/py/fastspecfit/emline_fit/sparse_rep.py +++ b/py/fastspecfit/emline_fit/sparse_rep.py @@ -1,8 +1,3 @@ -# -# Sparse representations of Jacobian of objective for emission line -# fitting (Ideal Jacobian rep is generated in EMLines_jacobian()) -# - import numpy as np from scipy.sparse.linalg import LinearOperator @@ -10,40 +5,54 @@ from .params_mapping import ParamsMapping -# -# Sparse Jacobian of objective function. For -# each camera's pixel range, the Jacobian -# is a matrix product -# -# W * M * J_I * J_S -# -# where -# J_S is the Jacobian of the parameter expansion -# J_I is the ideal Jacobian -# M is the camera's resolution matrix -# w is a diagonal matrix of per-observation weights -# -# Note that we precompute W*M*J_I for each camera -# with the external mulWMJ() function. This product -# has one contiguous run of nonzero entries per column, -# while J_S has either one or two nonzero entries -# per row. -# + class EMLineJacobian(LinearOperator): - - # - # CONSTRUCTOR ARGS: - # shape of Jacobian - # number of *line*-related free parameters in set - # array of start and end obs bin indices - # for each camera - # partial Jacobian jac = (W * M J_I) - # for each camera - # parameter expansion Jacobian J_S - # + """Sparse Jacobian of objective function for emission-line + fitting. For each camera's pixel range, the Jacobian is a matrix + product + + W * M * J_I * J_S + + where + J_S is the Jacobian of the parameter expansion + J_I is the ideal Jacobian for the Gaussian mixture + M is the camera's resolution matrix + W is a diagonal matrix of per-observation weights + + Note that we precompute W*M*J_I for each camera + with the external mulWMJ() function. This product + has one contiguous run of nonzero entries per column, + while J_S has either one or two nonzero entries + per row. + + When we use patch pedestals, the above Jacobian is + extended with additional rows corresponding to the + two free parameters (slope, intercept) of each patch. + + The components of this Jacobian are computed in params_mapping.py + for J_S and jacobian.py for everything else. + + """ + def __init__(self, shape, nLineFreeParms, camerapix, jacs, J_S, - J_P = None): + J_P = None): + """ + Parameters + ---------- + shape : :class:`tuple` [2] + Shape of Jacobian matrix. + nLineFreeParms : :class:`int` + Number of *line*-related free parameters in set. + camerapix : :class:`np.ndarray` [# cameras] + Array of start and end wavelength bin indices per camera. + jacs : :class:`tuple` [# cameras] + Tuple of partial Jacobians (W * M * J_I) for each camera. + J_S : :class:`np.ndarray` [# line params x # nLineFreeParms] + Parameter expansion Jacobian. + J_P : :class:`np.ndarray` [# patch params x # wavelength bins] or :None: + Jacobian of patch pedestals, if used. + """ self.camerapix = camerapix self.jacs = tuple(jacs) self.J_S = J_S @@ -60,11 +69,19 @@ def __init__(self, shape, nLineFreeParms, camerapix, jacs, J_S, super().__init__(dtype, shape) - # - # Compute left matrix-vector product J * v - # |v| = number of free parameters - # def _matvec(self, v): + """ + Compute left matrix-vector product J * v, where we are J. + + Parameters + ---------- + v : :class:`np.ndarray` [# of free parameters] + + Returns + ------- + w : :class:`np.ndarray [# of wavelength bins] + + """ nBins = self.shape[0] w = np.empty(nBins, dtype=v.dtype) @@ -88,13 +105,23 @@ def _matvec(self, v): return w - # - # Compute left matrix-matrix product J * M - # We do this just to avoid needing to ravel - # v in the more commonly used _matvec() - # def _matmat(self, M): + """ + Compute left matrix-matrix product J * M, where we are J. + This exists mainly to avoid having to ravel v in the + more frequently called _matvec(). M has an arbitrary + second dimension N. + + Parameters + ---------- + M : :class:`np.ndarray` [# of free parameters x N] + Returns + ------- + w : :class:`np.ndarray [N x # of wavelength bins] + + """ + nBins = self.shape[0] nVecs = M.shape[1] R = np.empty((nVecs, nBins), dtype=M.dtype) # transpose of result @@ -122,11 +149,19 @@ def _matmat(self, M): return R.T - # - # Compute right matrix product product v * J^T - # |v| = number of observable bins - # def _rmatvec(self, v): + """ + Compute right matrix-vector product v * J.T, where we are J. + + Parameters + ---------- + v : :class:`np.ndarray` [# of wavelength bins] + + Returns + ------- + w : :class:`np.ndarray [# of free parameters] + + """ nFreeParms = self.shape[1] w = np.zeros(nFreeParms, dtype=v.dtype) @@ -146,48 +181,48 @@ def _rmatvec(self, v): self._rmatvec_J(self.J_P, v, wPatches) return w - - # - # col_norms() - # Return the column norms of the Jacobian - # - def col_norms(self): - - norms = np.empty(self.shape[1], dtype=self.jacs[0][1].dtype) - - v = np.zeros(self.shape[0]) - for i in range(self.shape[1]): - - v[i] = 1. - - w = self._matvec(v) - norms[i] = np.norm(w) - - v[i] = 0. - - return norms # # Multiply ideal Jacobian J * v, writing result to w. # @staticmethod - @jit(nopython=True, fastmath=False, nogil=True) + @jit(nopython=True, nogil=True) def _matvec_J(J, v, w): + """ + Multiply partial Jacobian J with v, + returning result in supplied w. + Parameters + ---------- + J : :class:`np.ndarray` [# wavelength bins x # line parameters] + v : :class:`np.ndarray` [# of line parameters] + w : :class:`np.ndarray` [# wavength bins] (output param) + + """ + nbins = len(w) for j in range(nbins): w[j] = 0. _matvec_J_add(J, v, w) - # - # Multiply ideal Jacobian v * J^T, writing result to w. - # + @staticmethod - @jit(nopython=True, fastmath=False, nogil=True) + @jit(nopython=True, nogil=True) def _rmatvec_J(J, v, w): - + """ + Multiply v with partial Jacobian J.T, + returning result in supplied w. + + Parameters + ---------- + J : :class:`np.ndarray` [# wavelength bins x # line parameters] + v : :class:`np.ndarray` [# wavength bins] + w : :class:`np.ndarray` [# of line parameters] (output param) + + """ + endpts, values = J nvars = endpts.shape[0] @@ -200,37 +235,27 @@ def _rmatvec_J(J, v, w): acc += vals[j] * v[j + s] w[i] = acc - - # - # for debugging -- compute singular values of Jacobian - # and estimate its condition number. - # - def estimateConditionNumber(self): - - from scipi.sparse.linalg import svds - - try: - svs = svds(self, return_singular_vectors=False, - k=nFreeParms - 1, which="LM") - sv0 = svds(self, return_singular_vectors=False, - k=1, which="SM")[0] - cond = svs[-1] / sv0 - print(np.hstack((sv0, svs))) - print(f"cond(J) = {cond:.3e}") - - except: - print("Failed to compute Jacobian condition number") - - -# -# Multiply ideal Jacobian J * v, *adding* result to w. We have to -# make this a standalone function, rather than a static method of the -# EMLine_Jacobian class, because Numba cannot handle calling a static -# class method from JIT'd code. -# -@jit(nopython=True, fastmath=False, nogil=True) + +@jit(nopython=True, nogil=True) def _matvec_J_add(J, v, w): + """ + Multiply partial Jacobian J with v, + *adding* result to supplied w. + + Parameters + ---------- + J : :class:`np.ndarray` [# wavelength bins x # line parameters] + v : :class:`np.ndarray` [# of line parameters] + w : :class:`np.ndarray` [# wavength bins] (output param) + + Notes + ----- + This function is standalone, rather than a static method of + the EMLineJacobian class, only because Numba cannot call a static + class method from JIT'd code. + """ + endpts, values = J nvars = endpts.shape[0] diff --git a/py/fastspecfit/qa.py b/py/fastspecfit/qa.py index 16c7b8de..afdfaaa1 100644 --- a/py/fastspecfit/qa.py +++ b/py/fastspecfit/qa.py @@ -429,7 +429,7 @@ def major_formatter(x, pos): desismoothcontinuum.append(fullsmoothcontinuum[campix[0]:campix[1]]) # full model spectrum - _desiemlines = EMFit.emlinemodel_bestfit(fastspec, fastspec['Z'], np.hstack(data['wave']), data['res_emline'], + _desiemlines = EMFit.emlinemodel_bestfit(fastspec, fastspec['Z'], np.hstack(data['wave']), data['res'], data['camerapix'], snrcut=emline_snrmin) desiemlines = [] for icam in range(len(data['cameras'])): @@ -888,7 +888,7 @@ def major_formatter(x, pos): parameters[EMFit.doublet_idx] *= parameters[EMFit.doublet_src] lineprofiles = EMLine_MultiLines(parameters, np.hstack(data['wave']), redshift, EMFit.line_table['restwave'].value, - data['res_emline'], data['camerapix']) + data['res'], data['camerapix']) for iax, (minwave, maxwave, meanwave, deltawave, sig, linename, _linename) in enumerate( zip(minwaves, maxwaves, meanwaves, deltawaves, sigmas, linenames, _linenames)): From 4a0a6edb69172f68b643e2ac26f488fa13688ec6 Mon Sep 17 00:00:00 2001 From: jdbuhler Date: Tue, 3 Sep 2024 10:58:50 -0400 Subject: [PATCH 07/20] Update docs for jacobian.py to current std --- py/fastspecfit/emline_fit/jacobian.py | 77 +++++++++++++++++++-------- 1 file changed, 56 insertions(+), 21 deletions(-) diff --git a/py/fastspecfit/emline_fit/jacobian.py b/py/fastspecfit/emline_fit/jacobian.py index 35469b1f..fc055480 100644 --- a/py/fastspecfit/emline_fit/jacobian.py +++ b/py/fastspecfit/emline_fit/jacobian.py @@ -11,25 +11,40 @@ norm_pdf ) -# -# emline_model_jacobian() -# -# Compute the Jacobian of the function computed in emlines_model(). -# Inputs are as for emlines_model(). -# -# RETURNS: -# Sparse Jacobian as tuple (endpts, dd), where -# column j has nonzero values in interval [ endpts[j,0] , endpts[j,1] ) -# which are stored in dd[j]. -# -@jit(nopython=True, fastmath=False, nogil=True) + +@jit(nopython=True, nogil=True) def emline_model_jacobian(line_parameters, log_obs_bin_edges, ibin_widths, redshift, line_wavelengths, padding): - + """ + Compute Jacobian of the function computed in emline_model(). + + Parameters + ---------- + line_parameters : :class:`np.ndarray` + Parameters of all fitted lines. + log_obs_bin_edges : :class:`np.ndarray` [# wavelength bins + 1] + Natural logs of observed wavelength bin edges. + ibin_widths : :class:np.ndarray` [# wavelength bins] + Inverse of width of each observed wavelength bin. + redshift : :class:`np.float64` + Red shift of observed spectrum. + line_wavelengths : :class:`np.ndarray` [# lines] + Array of nominal wavelengths for all fitted lines. + padding : :class:`int`: + Number of entries to be added to size of each sparse row + allocated for output Jacobian, for later use. + + Returns + ------- + :class:`tuple` (endpts, dd) + Sparse Jacobian, where column j has nonzero values in interval + [ endpts[j,0] , endpts[j,1] ], which are stored in dd[j]. + + """ SQRT_2PI = np.sqrt(2*np.pi) @@ -165,19 +180,38 @@ def emline_model_jacobian(line_parameters, return (endpts, dd) -# compute partial Jacobian associated with just the patch parameters -# in the sparse column form used in sparse_rep.py. This Jacobian -# is independent of both the line model and the particular choices -# of slope/intercept for each patch, so can be computed once when -# we set up the optimization. -# -# FIXME: should we do this for each camera and apply resolution? + @staticmethod -@jit(nopython=True, fastmath=False, nogil=True) +@jit(nopython=True, nogil=True) def patch_jacobian(obs_bin_centers, obs_weights, patch_endpts, patch_pivotwave): + """ + Compute partial Jacobian associated with just the patch parameters + in the sparse column form used in sparse_rep.py. This Jacobian + is independent of both the line model and the particular choices + of slope/intercept for each patch, so can be computed once when + we set up the optimization. + + Parameters + ---------- + obs_bin_centers : :class:`np.ndarray` [# wavelength bins] + Center of each observed wavelength bin. + obs_weights : :class:`np.ndarray` [# wavelength bins] + Weights for each observed wavelength bin. + patch_endpts : :class:`np.ndarray` [# patches x 2] + Endpoints of each patch in wavelength array. + patch_pivotwave : :class:`np.ndarray` [# patches] + Wavelength offset for fitted affine params of each patch . + + Returns + ------- + :class:`tuple` (endpts, M) + Sparse Jacobian, where column j has nonzero values in interval + [ endpts[j,0] , endpts[j,1] ], which are stored in M[j]. + + """ nPatches = patch_endpts.shape[0] @@ -204,6 +238,7 @@ def patch_jacobian(obs_bin_centers, # patch. These derivatives are nonzero only # within the boundaries of the patch. # + M = np.empty((2*nPatches, maxPatchWidth)) for i in range(nPatches): s, e = endpts[i] From 516aafaea24e3e2e8b96f207bc6d8da9d25b742b Mon Sep 17 00:00:00 2001 From: jdbuhler Date: Tue, 3 Sep 2024 14:22:33 -0400 Subject: [PATCH 08/20] more docstring updates --- py/fastspecfit/emline_fit/model.py | 162 ++++++++++------- py/fastspecfit/emline_fit/params_mapping.py | 184 ++++++++++++++------ 2 files changed, 233 insertions(+), 113 deletions(-) diff --git a/py/fastspecfit/emline_fit/model.py b/py/fastspecfit/emline_fit/model.py index 079709ac..c7ba4611 100644 --- a/py/fastspecfit/emline_fit/model.py +++ b/py/fastspecfit/emline_fit/model.py @@ -10,31 +10,37 @@ norm_cdf ) -# -# emline_model() -# -# Given a fixed set of spectral lines and known redshift, and estimates -# for the amplitude, width, and velocity shift of each line, compute the -# average flux values observed in each of a series of spectral bins -# whose log wavelengths are bounded by log_obs_bin_edges. -# -# INPUTS: -# line_wavelengths -- array of nominal wavelengths for all fitted lines -# line_parameters -- parameters of each fitted line -# -# log_obs_bin_edges -- natural logs of observed wavelength bin edges -# redshift -- red shift of observed spectrum -# ibin_widths -- one over widths of each observed wavelength bin -# -# RETURNS: -# vector of average fluxes in each observed wavelength bin -# -@jit(nopython=True, fastmath=False, nogil=True) +@jit(nopython=True, nogil=True) def emline_model(line_wavelengths, line_parameters, log_obs_bin_edges, redshift, ibin_widths): + """ + Given a fixed set of spectral lines and known redshift, and estimates + for the amplitude, width, and velocity shift of each line, compute the + average flux values observed in each of a series of spectral bins + whose log wavelengths are bounded by log_obs_bin_edges. + + Parameters + ---------- + line_wavelengths : :class:`np.ndarray` [# lines] + Array of nominal wavelengths for all fitted lines. + line_parameters : :class:`np.ndarray` + Parameters of each fitted line. + log_obs_bin_edges : :class:`np.ndarray` [# obs wavelength bins + 1] + Natural logs of observed wavelength bin edges + redshift : :class:`np.float64` + Red shift of observed spectrum + ibin_widths : :class:`np.ndarray` [# obs wavelength bins] + Inverse widths of each observed wavelength bin. + + Returns + ------- + `np.ndarray` of length [# obs wavelength bins] with + average fluxes in each observed wavelength bin + + """ line_amplitudes, line_vshifts, line_sigmas = \ np.split(line_parameters, 3) @@ -70,23 +76,46 @@ def emline_model(line_wavelengths, return model_fluxes -# -# emline_perline_models() -# Given parameters for a set of lines, compute for each -# line individually its waveform in the range of bins -# described by log_bin_edges. The waveforms are returned -# sparsely in the same forma as the Jacobian below. -# -# This function shares the core computation of emline_model() -# but does not collapse the lines into one composite. -# -@jit(nopython=True, fastmath=False, nogil=True) +@jit(nopython=True, nogil=True) def emline_perline_models(line_wavelengths, line_parameters, log_obs_bin_edges, redshift, ibin_widths, padding): + """ + Given parameters for a set of lines, compute for each + line individually its waveform in the range of bins + described by log_bin_edges. The waveforms are returned + sparsely in the same forma as the Jacobian below. + + This function shares the core computation of emline_model() + but does not collapse the lines into one composite. + + Parameters + ---------- + line_wavelengths : :class:`np.ndarray` [# lines] + Array of nominal wavelengths for all fitted lines. + line_parameters : :class:`np.ndarray` + Parameters of each fitted line. + log_obs_bin_edges : :class:`np.ndarray` [# obs wavelength bins + 1] + Natural logs of observed wavelength bin edges + redshift : :class:`np.float64` + Red shift of observed spectrum + ibin_widths : :class:`np.ndarray` [# obs wavelength bins] + Inverse widths of each observed wavelength bin. + padding : :class:`int` + Number of entries by which to pad each sparse row for later use. + + Returns + ------- + Row-sparse matrix of size [# lines x # obs wavelength bins] + as tuple (endpts, profiles). Endpts[j] gives the + two endpoints of the range of wavelength bins with nonzero flux + for line j, while profiles[j] gives these nonzero values + (and may have some empty cells thereafter if padding is non-zero). + + """ line_amplitudes, line_vshifts, line_sigmas = \ np.split(line_parameters, 3) @@ -127,37 +156,7 @@ def emline_perline_models(line_wavelengths, return (endpts, line_profiles) -# -# emline_model_core() -# Compute the flux contribution of one spectral line to a model -# spectrum. We compute the range of bins where the line -# contributes flux, then write those contributions to the -# output array vals and return half-open range of bins -# [s, e) to which it contributes. -# -# INPUTS: -# line_wavelength -- nominal wavelength of line -# line_amplitude -- amplitude for spectral line -# line_vshift -- additional velocity shift for line -# line_sigma -- width of Gaussian profile for line -# -# log_obs_bin_edges -- natural log of observed wavelength bin edges -# redshift -- red shift of observed spectrum -# ibin_widths -- one over widths of each observed wavelength bin -# -# vals -- array in which to record values of -# any nonzero bin fluxes for this line -# -# RETURNS: -# tuple (s,e) such that nonzero fluxes occur in bins [s,e) -# (if s == e, no nonzero fluxes were computed) -# vals[0:e-s] contains these fluxes -# -# NB: vals MUST have length at least two more than the largest number -# of bins that could be computed. Entries in vals after the returned -# range may be set to arbitrary values. -# -@jit(nopython=True, fastmath=False, nogil=True) +@jit(nopython=True, nogil=True) def emline_model_core(line_wavelength, line_amplitude, line_vshift, @@ -166,6 +165,45 @@ def emline_model_core(line_wavelength, redshift, ibin_widths, vals): + """ + Compute the flux contribution of one spectral line to a model + spectrum. We compute the range of bins where the line + contributes flux, then write those contributions to the + output array vals and return half-open range of bins + [s, e) to which it contributes. + + Parameters + ---------- + line_wavelength : :class:`np.float64` + Nominal wavelength of line. + line_amplitude : :class:`np.float64` + Amplitude of line. + line_vshift : :class:`np.float64` + Velocity shift of line. + line_sigma : :class:`np.float64` + Width of Gaussian profile for line + log_obs_bin_edges : :class:`np.ndarray` [# obs wavelength bins + 1] + Natural logs of observed wavelength bin edges + redshift : :class:`np.float64` + Red shift of observed spectrum + ibin_widths : :class:`np.ndarray` [# obs wavelength bins] + Inverse widths of each observed wavelength bin. + vals : :class: `np.ndarray` (output) + Array in which to write values of any nonzero bin fluxes for this line + + Returns + ------- + Tuple (s,e) such that all nonzero fluxxes occur in the half-open + range of bins [s,e). If s == e, no nonzero fluxes were found. + The output parameter vals[0:e-s] contains these fluxes. + + Notes + ----- + vals MUST have length at least two more than the largest number + of bins that could be computed. Entries in vals after the returned + range may be set to arbitrary values. + + """ SQRT_2PI = np.sqrt(2 * np.pi) diff --git a/py/fastspecfit/emline_fit/params_mapping.py b/py/fastspecfit/emline_fit/params_mapping.py index 561d13d4..e299bf0c 100644 --- a/py/fastspecfit/emline_fit/params_mapping.py +++ b/py/fastspecfit/emline_fit/params_mapping.py @@ -1,20 +1,39 @@ -# -# Compute a mapping from the free parameters of a spectrum fitting problem -# to the full set of parameters for the EMLine model, as well the Jacobian -# of this mapping. We capture most of the complexity of the mapping once and -# do the minimum necessary updating to it for each new set of freeParameters. -# - import numpy as np from numba import jit class ParamsMapping(object): + """Compute a mapping from the free (line) parameters of a spectrum + fitting problem to the full set of parameters for the EMLine + model, as well the Jacobian of this mapping. We capture most of + the complexity of the mapping once and do the minimum necessary + work to update it for each new set of free parameters. + """ + def __init__(self, nParms, isFree, tiedSources, tiedFactors, doubletTargets, doubletSources): + """ + Parameters + ---------- + nParms : :class:`int` + Total number of parameters (free + fixed/tied). + isFree : :class:`np.ndarray` of `bool` [nParms] + Is each parameter free or fixed/tied? + tiedSources : :class:`np.ndarray` of `int` [nParms] + For parameters tied to a source param, idx of source + (-1 if not tied). + tiedFactors : :class:`np.ndarray` of `np.float64` [nParms] + For parameters tied to a source param, multiplier from + source to target. + doubletTargets: :class:`np.ndarray` of `int` + Parameters that are doublet ratios vs. a source. + doubletSources: :class:`np.ndarray` of `int` + Source parameters for all doublet ratios. + + """ self.nParms = nParms self.nFreeParms = np.sum(isFree) @@ -31,19 +50,37 @@ def __init__(self, nParms, self._precomputeJacobian() - - # return Boolean mask with "True" for each fixed parameter def fixedMask(self): + """ + Return Boolean mask with "True" for each fixed parameter. + """ return self.isFixed - # - # mapFreeToFull() - # Given a vector of free parameters, return the corresponding - # list of full parameters, accounting for fixed, tied, and - # doublet features. - # def mapFreeToFull(self, freeParms, out=None, patchDoublets=True): + """ + Given a vector of free parameters, return the corresponding + list of full parameters, accounting for fixed, tied, and + (if patchDoublets is true) doublet features. + + Parameters + ---------- + freeParms : :class:`np.ndarray` of `np.float64` + Values of free parameters. + out : :class: `np.ndarray` of `np.float64` or :None: + Vector to hold values of full parameter set; if none, + allocate one. + patchDoublets : :class: `bool` + Replace doublet target ratios by their actual values? + + Returns + ------- + A `np.ndarray` with the values of all parameters. + + Notes + ---- + Accelerated via Numba implementation below. + """ return self._mapFreeToFull(freeParms, self.nParms, @@ -53,14 +90,9 @@ def mapFreeToFull(self, freeParms, out=None, patchDoublets=True): out, patchDoublets) - # - # _mapFreeToFull() - # Given a vector of free parameters, return the corresponding - # list of full parameters, accounting for fixed, tied, and - # (if patchDoublets is true) doublet features. - # + @staticmethod - @jit(nopython=True, fastmath=False, nogil=True) + @jit(nopython=True, nogil=True) def _mapFreeToFull(freeParms, nParms, sources, factors, doubletPatches, fullParms, patchDoublets): @@ -78,27 +110,46 @@ def _mapFreeToFull(freeParms, nParms, sources, factors, return fullParms - # - # getJacobian() - # Given a vector v of free parameters, return the Jacobian - # of the transformation from free to full at v. The Jacobian - # is a sparse matrix represented as an array of nonzero entries - # (jacElts) and their values (jacFactors) - # def getJacobian(self, freeParms): + """ + Given a vector v of free parameters, return the Jacobian + of the transformation from free to full parameters at v. + + Parameters + ---------- + freeParms : :class:`np.ndarray` of `np.float64` + Vector of free parameter values + + Returns + ------- + Sparse Jacobian as a tuple (shape, elts, factors), where + - shape is shape of Jacobian [# parms x # free parms] + - elts gives the 2D indices of nonzero elements of Jacobian + - factors gives the values of the nonzero elements + + """ + for j, src_j_free in self.jacDoubletPatches: self.jacFactors[j] = freeParms[src_j_free] return ((self.nParms, self.nFreeParms), self.jacElts, self.jacFactors) - # - # Multiply parameter Jacobian J_S * v, writing result to w. - # @staticmethod - @jit(nopython=True, fastmath=False, nogil=True) + @jit(nopython=True, nogil=True) def _matvec(J_S, v, w): + """ + Multiply sparse parameter Jacobian J_S * v, writing result to w. + + Parameters + ---------- + J_S : :class: `tuple`: + Sparse Jacobian computed by getJacobian() + v : :class: `np.ndarray` [# free parameters] + w : :class: `np.ndarray` [# total parameters] (output) + """ + shape, elts, factors = J_S for j in range(shape[0]): # total params @@ -107,14 +158,22 @@ def _matvec(J_S, v, w): for i, (dst, src) in enumerate(elts): w[dst] += factors[i] * v[src] - - - # - # Multiply parameter Jacobian v * J_S^T, *adding* result to w. - # + @staticmethod - @jit(nopython=True, fastmath=False, nogil=True) + @jit(nopython=True, nogil=True) def _add_rmatvec(J_S, v, w): + """ + Multiply v with sparse parameter Jacobian J_S.T, + *adding* result to w. + + Parameters + ---------- + J_S : :class: `tuple`: + Sparse Jacobian computed by getJacobian() + v : :class: `np.ndarray` [# total parameters] + w : :class: `np.ndarray` [# free parameters] (output) + + """ _, elts, factors = J_S @@ -124,14 +183,34 @@ def _add_rmatvec(J_S, v, w): ########################################################### - # - # Precompute all the transformations from free parameters - # to full parameters that do not require knowledge of the - # free parameter values. - # def _precomputeMapping(self, isFree, pFree, tiedSources, tiedFactors, doubletTargets, doubletSources): + """ + Precompute and store all the transformations from free + parameters to full parameters that do not require knowledge of + the free parameter values. Create "patches" for doublets so + that we can quickly compute the values of all doublet targets + once their source values and ratios are known. + + isFree : :class:`np.ndarray` of `bool` [nParms] + Is each parameter free or fixed/tied? + pFree: :class:`np.ndarray` of `int` [nParms] + If a parameter in the full list isFree, which + index in the free parameter list does it come from? + (If not free, entry is undefined) + tiedSources : :class:`np.ndarray` of `int` [nParms] + For parameters tied to a source param, idx of source + (-1 if not tied). + tiedFactors : :class:`np.ndarray` of `np.float64` [nParms] + For parameters tied to a source param, multiplier from + source to target. + doubletTargets: :class:`np.ndarray` of `int` + Parameters that are doublet ratios vs. a source. + doubletSources: :class:`np.ndarray` of `int` + Source parameters for all doublet ratios. + + """ # by default, assume parameters are fixed sources = np.full(self.nParms, -1, dtype=np.int32) @@ -170,14 +249,17 @@ def _precomputeMapping(self, isFree, pFree, self.isFixed = (self.sources == -1) - # - # Precompute as much of the Jacobian of the transformation - # from free parameters to full parameters as does not require - # knowledge of the free parameter values. We rely on the - # precomputation for the mapping to avoid recalculating all - # the source/factor information. - # def _precomputeJacobian(self): + """ + Precompute and store as much of the Jacobian of the + transformation from free parameters to full parameters as does + not require knowledge of the free parameter values. We rely + on the precomputation for the mapping to avoid recalculating + all the source/factor information. We transform doublet + patches for the mapping into a larger set of patches for the + Jacobian. + + """ isLive = np.logical_not(self.isFixed) liveIdx = np.where(isLive)[0] From 2df2666829d32f95c9bf2986fdc0b950c4200b20 Mon Sep 17 00:00:00 2001 From: jdbuhler Date: Tue, 3 Sep 2024 14:25:19 -0400 Subject: [PATCH 09/20] Reduce least_squares overhead by switching from lsmr to exact TRF, and further limit ftol to 1e-3 for greater speed. --- doc/changes.rst | 4 ++++ py/fastspecfit/continuum.py | 4 ++-- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/doc/changes.rst b/doc/changes.rst index fd89c862..eafd105b 100644 --- a/doc/changes.rst +++ b/doc/changes.rst @@ -5,6 +5,9 @@ Change Log 3.0.0 (not released yet) ------------------------ +* Use Numba-based Resolution implementation with faster multiply. + Also tweak parameters of continuum least_squares optimizer for + speed. [`PR #181`_] * Cache Fourier transforms of the templates for speed [`PR #180`_]. * Update line-list (drop broad HeI lines) and rest wavelengths [`PR #179`_]. * Near-total rewrite of both the continuum and emission-line fitting engines and @@ -14,6 +17,7 @@ Change Log .. _`PR #177`: https://github.com/desihub/fastspecfit/pull/177 .. _`PR #179`: https://github.com/desihub/fastspecfit/pull/179 .. _`PR #180`: https://github.com/desihub/fastspecfit/pull/180 +`` _`PR #181`: https://github.com/desihub/fastspecfit/pull/181 2.5.2 (2024-04-28) ------------------ diff --git a/py/fastspecfit/continuum.py b/py/fastspecfit/continuum.py index b7b9abc0..8d18db9e 100644 --- a/py/fastspecfit/continuum.py +++ b/py/fastspecfit/continuum.py @@ -983,8 +983,8 @@ def fit_stellar_continuum(self, templateflux, fit_vdisp, conv_pre=None, # convergence. There may be faster ways, of course... fit_info = least_squares(self._stellar_objective, initial_guesses, kwargs=farg, bounds=tuple(zip(*bounds)), method='trf', - tr_solver='lsmr', tr_options={'regularize': True}, - x_scale='jac', max_nfev=5000, ftol=1e-5, xtol=1e-10)#, verbose=2) + tr_solver='exact', tr_options={'regularize': True}, + x_scale='jac', max_nfev=5000, ftol=1e-3, xtol=1e-10)#, verbose=2) bestparams = fit_info.x resid = fit_info.fun From d5d8f64185b3b0127ace856d78bcc0d76989199c Mon Sep 17 00:00:00 2001 From: jdbuhler Date: Tue, 3 Sep 2024 14:44:10 -0400 Subject: [PATCH 10/20] * fix markup error in changes.rst * attempt to clean up some formatting warnings * add new-style docstrings to emline_fit/utils.py --- doc/changes.rst | 2 +- py/fastspecfit/emline_fit/jacobian.py | 13 +++-- py/fastspecfit/emline_fit/model.py | 3 +- py/fastspecfit/emline_fit/params_mapping.py | 14 +++--- py/fastspecfit/emline_fit/sparse_rep.py | 8 +-- py/fastspecfit/emline_fit/utils.py | 56 ++++++++++++--------- 6 files changed, 54 insertions(+), 42 deletions(-) diff --git a/doc/changes.rst b/doc/changes.rst index eafd105b..fceba9dc 100644 --- a/doc/changes.rst +++ b/doc/changes.rst @@ -17,7 +17,7 @@ Change Log .. _`PR #177`: https://github.com/desihub/fastspecfit/pull/177 .. _`PR #179`: https://github.com/desihub/fastspecfit/pull/179 .. _`PR #180`: https://github.com/desihub/fastspecfit/pull/180 -`` _`PR #181`: https://github.com/desihub/fastspecfit/pull/181 +.. _`PR #181`: https://github.com/desihub/fastspecfit/pull/181 2.5.2 (2024-04-28) ------------------ diff --git a/py/fastspecfit/emline_fit/jacobian.py b/py/fastspecfit/emline_fit/jacobian.py index fc055480..55e289a9 100644 --- a/py/fastspecfit/emline_fit/jacobian.py +++ b/py/fastspecfit/emline_fit/jacobian.py @@ -49,21 +49,21 @@ def emline_model_jacobian(line_parameters, SQRT_2PI = np.sqrt(2*np.pi) nbins = len(log_obs_bin_edges) - 1 - + line_amplitudes, line_vshifts, line_sigmas = \ np.split(line_parameters, 3) - + # buffers for per-parameter calculations, sized large # enough for max possible range [s .. e), plus extra # padding as directed by caller. max_width = max_buffer_width(log_obs_bin_edges, line_sigmas, padding) - + nlines = len(line_wavelengths) dd = np.empty((3 * nlines, max_width), dtype=line_amplitudes.dtype) endpts = np.zeros((3 * nlines, 2), dtype=np.int32) - starts = endpts[:,0] - ends = endpts[:,1] + starts = endpts[:, 0] + ends = endpts[:, 1] # compute partial derivatives for avg values of all Gaussians # inside each bin. For each Gaussian, we only compute @@ -94,7 +94,7 @@ def emline_model_jacobian(line_parameters, if hi == 0 or lo == len(log_obs_bin_edges): continue - nedges = hi - lo + 2 # compute values at edges [lo - 1 ... hi] + nedges = hi - lo + 2 # compute values at edges [lo - 1 ... hi] # Compute contribs of each line to each partial derivative in place. # No sharing of params between peaks means that we never have to @@ -180,7 +180,6 @@ def emline_model_jacobian(line_parameters, return (endpts, dd) - @staticmethod @jit(nopython=True, nogil=True) def patch_jacobian(obs_bin_centers, diff --git a/py/fastspecfit/emline_fit/model.py b/py/fastspecfit/emline_fit/model.py index c7ba4611..baf2006a 100644 --- a/py/fastspecfit/emline_fit/model.py +++ b/py/fastspecfit/emline_fit/model.py @@ -10,6 +10,7 @@ norm_cdf ) + @jit(nopython=True, nogil=True) def emline_model(line_wavelengths, line_parameters, @@ -239,7 +240,7 @@ def emline_model_core(line_wavelength, # vals[i] --> edge i + lo - 1 - vals[0] = 0. # edge lo - 1 + vals[0] = 0. # edge lo - 1 for i in range(1, nedges-1): diff --git a/py/fastspecfit/emline_fit/params_mapping.py b/py/fastspecfit/emline_fit/params_mapping.py index e299bf0c..527ef15c 100644 --- a/py/fastspecfit/emline_fit/params_mapping.py +++ b/py/fastspecfit/emline_fit/params_mapping.py @@ -2,6 +2,7 @@ from numba import jit + class ParamsMapping(object): """Compute a mapping from the free (line) parameters of a spectrum fitting problem to the full set of parameters for the EMLine @@ -99,11 +100,11 @@ def _mapFreeToFull(freeParms, nParms, sources, factors, for j, src_j_free in doubletPatches: factors[j] = freeParms[src_j_free] if patchDoublets else 1. - if fullParms == None: + if fullParms is None: fullParms = np.empty(nParms, dtype=freeParms.dtype) for j, src_j_free in enumerate(sources): - fullParms[j] = factors[j] # copy fixed zeros + fullParms[j] = factors[j] # copy fixed zeros if src_j_free != -1: fullParms[j] *= freeParms[src_j_free] @@ -152,7 +153,7 @@ def _matvec(J_S, v, w): shape, elts, factors = J_S - for j in range(shape[0]): # total params + for j in range(shape[0]): # total params w[j] = 0. for i, (dst, src) in enumerate(elts): @@ -180,9 +181,10 @@ def _add_rmatvec(J_S, v, w): for i, (dst, src) in enumerate(elts): w[src] += factors[i] * v[dst] - + ########################################################### - + + def _precomputeMapping(self, isFree, pFree, tiedSources, tiedFactors, doubletTargets, doubletSources): @@ -293,7 +295,7 @@ def _precomputeJacobian(self): # jacElts already has coefficient (j, p[j]). # its factor should be patched from v[ p[src_j] ] - jacDoubletPatches[2*i, :] = (liveOffsets[j], p_src_j) + jacDoubletPatches[2*i, :] = (liveOffsets[j], p_src_j) # add a second coefficient (j, p[src_j]). # its factor should be patched from v[ p[j] ] diff --git a/py/fastspecfit/emline_fit/sparse_rep.py b/py/fastspecfit/emline_fit/sparse_rep.py index 02a8ec3b..55ec92ba 100644 --- a/py/fastspecfit/emline_fit/sparse_rep.py +++ b/py/fastspecfit/emline_fit/sparse_rep.py @@ -34,8 +34,8 @@ class EMLineJacobian(LinearOperator): """ - def __init__(self, shape, nLineFreeParms, camerapix, jacs, J_S, - J_P = None): + def __init__(self, shape, nLineFreeParms, camerapix, jacs, + J_S, J_P=None): """ Parameters ---------- @@ -63,7 +63,7 @@ def __init__(self, shape, nLineFreeParms, camerapix, jacs, J_S, nPatchParms = J_P[0].shape[0] dtype = jacs[0][1].dtype - nParms = jacs[0][1].shape[0] # num line params in full set + nParms = jacs[0][1].shape[0] # num line params in full set self.vFull = np.empty(nParms, dtype=dtype) super().__init__(dtype, shape) @@ -124,7 +124,7 @@ def _matmat(self, M): nBins = self.shape[0] nVecs = M.shape[1] - R = np.empty((nVecs, nBins), dtype=M.dtype) # transpose of result + R = np.empty((nVecs, nBins), dtype=M.dtype) # transpose of result for i in range(nVecs): w = R[i,:] diff --git a/py/fastspecfit/emline_fit/utils.py b/py/fastspecfit/emline_fit/utils.py index cac272ed..2813de47 100644 --- a/py/fastspecfit/emline_fit/utils.py +++ b/py/fastspecfit/emline_fit/utils.py @@ -5,36 +5,37 @@ from fastspecfit.util import C_LIGHT + # Do not bother computing normal PDF/CDF if more than this many # standard deviations from mean. MAX_SDEV = 5. -# -# norm_pdf() -# PDF of standard normal distribution at a point a -# -@jit(nopython=True, fastmath=False, nogil=True) -def norm_pdf(a): +@jit(nopython=True, nogil=True) +def norm_pdf(a): + """ + PDF of standard normal distribution at a point a + """ + SQRT_2PI = np.sqrt(2 * np.pi) return 1/SQRT_2PI * np.exp(-0.5 * a**2) -# -# norm_cdf() -# Approximate the integral of a standard normal PDF from -infty to a. -# -# Optimization (currently disabled because it is not needed): If -# |a| > MAX_SDEV, treat the value as extreme and return 0 or 1 as -# appropriate. -# -@jit(nopython=True, fastmath=False, nogil=True) + +@jit(nopython=True, nogil=True) def norm_cdf(a): + """ + Approximate the integral of a standard normal PDF from -infty to a. + """ SQRT1_2 = 1.0 / np.sqrt(2) z = np.abs(a) + # Optimization (currently disabled because it is not needed): If + # |a| > MAX_SDEV, treat the value as extreme and return 0 or 1 as + # appropriate. + #if z > MAX_SDEV: # short-circuit extreme values # if a > 0: # y = 1. @@ -50,16 +51,25 @@ def norm_cdf(a): return y -# -# max_buffer_width() -# Compute a safe estimate of the number of nonzero bin fluxes possible -# for a line spanning a subrange of bins with edges log_obs_bin_edges, -# assuming the line's width is one of the values in line_sigmas. -# Optionally add 2*padding to allow future expansion to left and right. -# -@jit(nopython=True, fastmath=False, nogil=True) +@jit(nopython=True, nogil=True) def max_buffer_width(log_obs_bin_edges, line_sigmas, padding=0): + """ + Compute a safe estimate of the number of nonzero bin fluxes possible + for a line spanning a subrange of bins with edges log_obs_bin_edges, + assuming the line's width is one of the values in line_sigmas. + Optionally add 2*padding to allow future expansion to left and right. + Parameters + ---------- + log_obs_bin_edges : :class:`np.ndarray` [# obs wavelength bins + 1] + log of wavelengths of all observed bin edges. + line_sigmas : :class:`np.ndarray` [# nlines] + Gaussian widths of all spectral lines. + padding : :class:`int` + Padding parameter to add to width for future use. + + """ + # Find the largest width sigma for any line, and # allocate enough space for twice that much width # in bins, given the smallest observed bin width. From 267571a2d70917cd04565f8bebd0431c9f415bf5 Mon Sep 17 00:00:00 2001 From: jdbuhler Date: Tue, 3 Sep 2024 16:16:42 -0400 Subject: [PATCH 11/20] update docstrings for interface, and make a couple of functions inner where they have only one caller --- py/fastspecfit/emline_fit/interface.py | 535 +++++++++++++++---------- 1 file changed, 332 insertions(+), 203 deletions(-) diff --git a/py/fastspecfit/emline_fit/interface.py b/py/fastspecfit/emline_fit/interface.py index 1a334dd4..8759d4b1 100644 --- a/py/fastspecfit/emline_fit/interface.py +++ b/py/fastspecfit/emline_fit/interface.py @@ -1,8 +1,3 @@ -# -# Calculation of objective function and -# Jacobian for emission line fitting -# - import numpy as np from math import erf, erfc @@ -30,7 +25,10 @@ class EMLine_Objective(object): - + """ + Objective function for emission-line fitting. + """ + def __init__(self, obs_bin_centers, obs_fluxes, @@ -41,6 +39,28 @@ def __init__(self, camerapix, params_mapping, continuum_patches=None): + """ + Parameters + ---------- + obs_bin_centers : :class:`np.ndarray` [# obs wavelength bins] + Center wavelength of each observed wavelength bin. + obs_fluxes : :class:`np.ndarray` [# obs wavelength bins] + Observed flux in each wavelength bin. + obs_weights : :class:`np.ndarray` [# obs wavelength bins] + Weighting factors for each wavelength bin in residual. + redshift : :class:`np.float64` + Red shift of observed spectrum. + line_wavelengths : :class:`np.ndarray` [# lines] + Array of nominal wavelengths for all fitted lines. + resolution_matrices : tuple of :class:`fastspecfit.resolution.Resolution` + Resolution matrices for each camera. + params_mapping : :class:`.params_mapping.ParamsMapping` + Mapping from free to full parameters. + continuum_patches : patch dictionary or :None: + If using patch pedestals, a dictionary of patch parameters; + else None. + + """ self.dtype = obs_fluxes.dtype @@ -74,15 +94,26 @@ def __init__(self, # on every call to objective/jacobian self.line_parameters = \ np.empty(self.params_mapping.nParms, dtype=self.dtype) - - # - # Objective function for least-squares optimization - # - # Build the emline model as described above and compute the weighted - # vector of residuals between the modeled fluxes and the observations. - # + + def objective(self, free_parameters): + """ + Objective function for least-squares optimization. + + Build the emline model as described above and compute the weighted + vector of residuals between the modeled fluxes and the observations. + Parameters + ---------- + free_parameters : :class:`np.ndarray` + Values of all free parameters. + + Returns + ------- + `np.ndarray` of residuals for each wavelength bin. + + """ + nLineParameters = len(free_parameters) - 2 * self.nPatches line_free_parameters = free_parameters[:nLineParameters] @@ -110,13 +141,12 @@ def objective(self, free_parameters): intercepts = free_parameters[nLineParameters+self.nPatches:] # add patch pedestals to line model - _add_patches(model_fluxes, - slopes, - intercepts, + _add_patches(self.obs_bin_centers, model_fluxes, self.patch_endpts, self.patch_pivotwave, - self.obs_bin_centers) - + slopes, + intercepts) + # turn model fluxes into residuals in-place to avoid # unwanted memory allocation residuals = model_fluxes @@ -126,14 +156,25 @@ def objective(self, free_parameters): return residuals - # - # Jacobian of objective function for least-squares EMLine - # optimization. The result of the detailed calculation is converted - # to a sparse matrix, since it is extremely sparse, to speed up - # subsequent matrix-vector multiplies in the optimizer. - # def jacobian(self, free_parameters): + """ + Jacobian of objective function for emission line fitting + optimization. The result of the detailed calculation is converted + to a sparse matrix, since it is extremely sparse, to speed up + subsequent matrix-vector multiplies in the optimizer. + Parameters + ---------- + free_parameters : :class:`np.ndarray` + Values of all free parameters. + + Returns + ------- + Sparse Jacobian at given parameter value as + defined in sparse_rep.py. + + """ + nLineFreeParms = len(free_parameters) - 2 * self.nPatches line_free_parameters = free_parameters[:nLineFreeParms] @@ -181,13 +222,11 @@ def jacobian(self, free_parameters): self.J_P) return J - + ################################################################################## -# -# compatibility entry point to compute modeled fluxes -# + def build_model(redshift, line_parameters, line_wavelengths, @@ -195,6 +234,28 @@ def build_model(redshift, resolution_matrices, camerapix, continuum_patches=None): + """ + Compatibility entry point to compute modeled fluxes. + + Parameters + ---------- + redshift : :class:`np.float64` + Red shift of observed spectrum. + line_parameters : :class:`np.ndarray` + Array of all fitted line parameters. + line_wavelengths : :class:`np.ndarray` [# lines] + Array of nominal wavelengths for all fitted lines. + obs_bin_centers : :class:`np.ndarray` [# obs wavelength bins] + Center wavelength of each observed wavelength bin. + resolution_matrices : tuple of :class:`fastspecfit.resolution.Resolution` + Resolution matrices for each camera. + camerapix : :class:`np.ndarray` of `int` [# cameras x 2] + Pixels corresponding to each camera in obs wavelength array. + continuum_patches : patch dictionary or :None: + If using patch pedestals, a dictionary of patch parameters; + else None. + + """ log_obs_bin_edges, ibin_widths = _prepare_bins(obs_bin_centers, camerapix) @@ -214,38 +275,23 @@ def build_model(redshift, if continuum_patches is not None: # add patch pedestals to model fluxes - _add_patches(model_fluxes, - continuum_patches['slope'].value, - continuum_patches['intercept'].value, + _add_patches(obs_bin_centers, model_fluxes, continuum_patches['endpts'].value, - continuum_patches['pivotwave'].value, - obs_bin_centers) + continuum_patches['pivotwave'].value, + continuum_patches['slope'].value, + continuum_patches['intercept'].value) return model_fluxes -# -# class MultiLines -# Construct models for each of a list of lines across -# one or more cameras. Given a line number, return -# a rendered model for that line as a sparse array. -# class MultiLines(object): - - # Given fitted parameters for all emission lines, - # compute models for each line for each camera - # - # INPUTS: - # parameters -- full array of fitted line parameters - # obs_bin_centers -- center wavelength of each - # observed flux bin - # redshift -- object redshift - # line_wavelengths -- wavelengths of all fitted lines - # resolution_matrices -- list of sparse resolution matrices - # in same form taken by objective - # camerapix -- wavelength ranges for each camera - # + """ + Construct models for each of a list of lines across one or + more cameras. Return an object that lets caller obtain + a rendered model for each individual line as a sparse array. + """ + def __init__(self, line_parameters, obs_bin_centers, @@ -253,7 +299,37 @@ def __init__(self, line_wavelengths, resolution_matrices, camerapix): - + """ + Given fitted parameters for all emission lines, compute + models for each line for each camera. + + Parameters + ---------- + line_parameters : :class:`np.ndarray` + Array of all fitted line parameters. + obs_bin_centers : :class:`np.ndarray` [# obs wavelength bins] + Center wavelength of each observed wavelength bin. + redshift : :class:`np.float64` + Red shift of observed spectrum. + line_wavelengths : :class:`np.ndarray` [# lines] + Array of nominal wavelengths for all fitted lines. + resolution_matrices : tuple of :class:`fastspecfit.resolution.Resolution` + Resolution matrices for each camera. + camerapix : :class:`np.ndarray` of `int` [# cameras x 2] + Pixels corresponding to each camera in obs wavelength array. + + """ + + @jit(nopython=True, nogil=True) + def _suppress_negative_fluxes(endpts, M): + """ + suppress negative fluxes arising from resolution matrix + """ + for i in range(M.shape[0]): + s, e = endpts[i] + for j in range(e-s): + M[i,j] = np.maximum(M[i,j], 0.) + self.line_models = [] _build_multimodel_core(line_parameters, obs_bin_centers, @@ -265,34 +341,30 @@ def __init__(self, self.line_models = tuple(self.line_models) - # suppress negative fluxes arising from resolution matrix for endpts, M in self.line_models: - self._suppress_negative_fluxes(endpts, M) - - @staticmethod - @jit(nopython=True, fastmath=False, nogil=True) - def _suppress_negative_fluxes(endpts, M): - for i in range(M.shape[0]): - s, e = endpts[i] - for j in range(e-s): - M[i,j] = np.maximum(M[i,j], 0.) - - # - # getLine(): - # Return a model for one emission line - # - # INPUTS: line -- number of line to return - # - # RETURNS: sparse line model as a tple - # (s, e), data - # - # where the range (s,e) in the combined bin array across - # all cameras (as in obs_bin_centers) contains all bins - # with nonzero fluxes for that line, and data is an array - # of length e - s that contains the bin values. - # + _suppress_negative_fluxes(endpts, M) + + + def getLine(self, line): + """ + Return a model for one emission line. + Parameters + ---------- + line : :class:`int` + Index of line to return. + + Returns + ------- + Sparse line model as tuple ((s, e), data), where the range + [s,e) in the combined bin array across all cameras (as in + obs_bin_centers) contains all bins with nonzero fluxes for + that line, and data is an array of length e - s that contains + the bin values. + + """ + s = 1000000000 e =-1 @@ -333,31 +405,56 @@ def getLine(self, line): return (s, e), data - -# find_peak_amplitudes() -# Given fitted parameters for all emission lines, report for -# each line the largest flux that it contributes to any -# observed bin. -# -# INPUTS: -# parameters -- full array of fitted line parameters -# bin_data -- preprocessed bin data in the same form -# taken by the optimizer objective -# redshift -- object redshift -# line_wavelengths -- wavelengths of all fitted lines -# resolution_matrices -- list of sparse resolution matrices -# in same form taken by objective -# camerapix -- wavelength ranges for each camera -# -# RETURNS: -# an array of maximum amplitudes for each line -# + def find_peak_amplitudes(line_parameters, obs_bin_centers, redshift, line_wavelengths, resolution_matrices, camerapix): + """ + Given fitted parameters for all emission lines, report for + each line the largest flux that it contributes to any observed + bin. + + Parameters + ---------- + parameters : :class:`np.ndarray` + Array of all fitted line parameters. + bin_data : :class:`tuple` + Preprocessed bin data in form returned by _prepare_bins + redshift : :class:`np.float64` + Red shift of observed spectrum. + line_wavelengths : :class:`np.ndarray` [# lines] + Array of nominal wavelengths for all fitted lines. + resolution_matrices : tuple of :class:`fastspecfit.resolution.Resolution` + Resolution matrices for each camera. + camerapix : :class:`np.ndarray` of `int` [# cameras x 2] + Pixels corresponding to each camera in obs wavelength array. + + Returns + ------- + Array with maximum amplitude observed for each line. + + """ + + @jit(nopython=True, nogil=True) + def _update_line_maxima(max_amps, line_models): + """ + Given an array of line waveforms and an array of maximum + amplitudes observed for each line so far, if the waveform contains + a higher amplitude than the previous maximum, update the maximum + in place. + """ + endpts, vals = line_models + + # find the highest flux for each peak; if it's + # bigger than any seen so far, update global max + for i in range(vals.shape[0]): + ps, pe = endpts[i] + if pe > ps: + max_amps[i] = np.maximum(max_amps[i], + np.max(vals[i,:pe-ps])) max_amps = np.zeros_like(line_wavelengths, dtype=line_parameters.dtype) @@ -371,52 +468,10 @@ def find_peak_amplitudes(line_parameters, return max_amps - -# -# update_line_maxima() -# Given an array of line waveforms and an array of -# maximum amplitudes observed for each line so far, -# if the waveform contains a higher amplitude than -# the previous maximum, update the maximum in place. -# -@jit(nopython=True, fastmath=False, nogil=True) -def _update_line_maxima(max_amps, line_models): - - endpts, vals = line_models - - # find the highest flux for each peak; if it's - # bigger than any seen so far, update global max - for i in range(vals.shape[0]): - ps, pe = endpts[i] - if pe > ps: - max_amps[i] = np.maximum(max_amps[i], - np.max(vals[i,:pe-ps])) - ########################################################################## -# -# _build_model_core() -# Core loop for computing a combined model flux -# from a set of spectral emission lines -# -# INPUTS: -# line_parameters - parameters of each line -# line_wavelengths - wavelength of each line -# redshift - object redshift -# log_obs_bin_edges - log-wavelengths of edges -# for each observed flux bin -# ibin_widths - 1/width of each flux bin -# rsolution_matrices - resolution matrices per camera -# camerapix - ranges of bins associated with -# each camera (2D nparray) -# model_fluxes - output array -# -# NB: log_obs_bin_edges and ibin_widths must be -# padded as is done by _prepare_bins() below -# -# RETURNS: combined flux model in model_fluxes -# + def _build_model_core(line_parameters, line_wavelengths, redshift, @@ -425,6 +480,30 @@ def _build_model_core(line_parameters, resolution_matrices, camerapix, model_fluxes): + """ + Core loop for computing a combined model flux from a set of + spectral emission lines. + + Parameters + ---------- + line_parameters : :class:`np.ndarray` + Combined array of amplitudes, vshifts, and sigmas + line_wavelengths : :class:`np.ndarray` [# lines] + Array of nominal wavelengths for all fitted lines. + redshift : :class:`np.float64` + Red shift of observed spectrum. + log_obs_bin_edges : :class:`np.ndarray` [# obs wavelength bins + 1] + Natural logs of observed wavelength bin edges + ibin_widths : :class:`np.ndarray` [# obs wavelength bins] + Inverse widths of each observed wavelength bin. + resolution_matrices : tuple of :class:`fastspecfit.resolution.Resolution` + Resolution matrices for each camera. + camerapix : :class:`np.ndarray` of `int` [# cameras x 2] + Pixels corresponding to each camera in obs wavelength array. + model_fluxes : :class:`np.ndarray` [# obs wavelength bins] (output) + Returns computed model flux for each wavelength bin. + + """ for icam, campix in enumerate(camerapix): @@ -448,23 +527,6 @@ def _build_model_core(line_parameters, resolution_matrices[icam].dot(mf, model_fluxes[s:e]) -# -# _build_multimodel_core() -# Core loop for computing individual line flux models -# sparsely from a set of spectral emission lines -# -# INPUTS: -# parameters - combined array of amplitudes, vshifts, and sigmas -# obs_bin_centers - center wavelength of each observed flux bin -# redshift - object redshift -# line_wavelengths - wavelength of each line -# rsolution_matrices - resolution matrices per camera -# camerapix - ranges of bins associated with -# each camera (2D nparray) -# consumer_fun - function to which we pass the -# computed array of line models -# for each camera -# def _build_multimodel_core(line_parameters, obs_bin_centers, redshift, @@ -472,6 +534,30 @@ def _build_multimodel_core(line_parameters, resolution_matrices, camerapix, consumer_fun): + """ + Core loop for computing array of individual line flux models sparsely + from a set of spectral emission lines. Result is not returned but + passed to a consumer function supplied by the caller. + + Parameters + ---------- + line_parameters : :class:`np.ndarray` + Combined array of amplitudes, vshifts, and sigmas + obs_bin_centers : :class:`np.ndarray` [# obs wavelength bins] + Center wavelength of each observed wavelength bin. + redshift : :class:`np.float64` + Redshift of observed spectrum. + line_wavelengths : :class:`np.ndarray` [# lines] + Array of nominal wavelengths for all fitted lines. + resolution_matrices : tuple of :class:`fastspecfit.resolution.Resolution` + Resolution matrices for each camera. + camerapix : :class:`np.ndarray` of `int` [# cameras x 2] + Pixels corresponding to each camera in obs wavelength array. + consumer_fun : + Function to which we pass computed sparse array of line models + for each camera. + + """ log_obs_bin_edges, ibin_widths = _prepare_bins(obs_bin_centers, camerapix) @@ -480,7 +566,6 @@ def _build_multimodel_core(line_parameters, # start and end for obs fluxes of camera icam s, e = campix - # Actual inverse bin widths are in ibin_widths[s+1:e+2]. # Setup guarantees that at least one more array entry @@ -507,19 +592,34 @@ def _build_multimodel_core(line_parameters, consumer_fun((endpts, M)) -# -# _add_patches() -# Add patch pedestals for patches with given -# endpts, slopes, intercepts, and pivots to -# an array of model fluxes -# -@jit(nopython=True, fastmath=False, nogil=True) -def _add_patches(model_fluxes, - slopes, - intercepts, +@jit(nopython=True, nogil=True) +def _add_patches(obs_bin_centers, + model_fluxes, patch_endpts, patch_pivotwaves, - obs_bin_centers): + slopes, + intercepts): + """ + Add patch pedestals for patches with given endpts, slopes, + intercepts, and pivots to an array of model fluxes. + + Parameters + ---------- + obs_bin_centers : :class:`np.ndarray` [# obs wavelength bins] + Center wavelength of each obesrved bin. + model_fluxes : :class:`np.ndarray` [# obs wavelength bins] (output) + Fluxes computed from base emission line model. Modified + by addition of patch pedestal contribution to each bin. + patch_endpts: :class:`np.ndarray` of `int` [2 * # patches] + Endpoint indices of each patch in wavelength bin array + patch_pivotwave: :class:`np.ndarray` [# patches] + Offset for each patch pedestal. + slopes : :class:`np.ndarray` [# patches] + Slope of each patch pedestal. + itercepts: :class:`np.ndarray` [# patches] + Intercept of each patch pedestal. + + """ # add patch pedestals to line model nPatches = len(slopes) @@ -538,24 +638,42 @@ def _add_patches(model_fluxes, ################################################################### -# -# mulWMJ() -# Compute the sparse matrix product P = WMJ, where -# W is a diagonal matrix (represented by a vector) -# M is a resolution matrix (the .data field of a -# ResMatrix, since Numba can't handle classes) -# J is a column-sparse matrix giving the nonzero -# values in one contiguous range per column -# -# We assume that J jac has enough free space in -# each column to store the results of the multiply -# and update it in place, returning the same -# jac as our return value. -# -@jit(nopython=True, fastmath=False, nogil=True) -def mulWMJ(w, M, jac): - - endpts, J = jac + +@jit(nopython=True, nogil=True) +def mulWMJ(w, M, Jsp): + """Compute the sparse matrix product P = WMJ, where + W is a diagonal weight matrix + M is a resolution matrix + J is a column-sparse matrix giving the nonzero + values in one contiguous range per column + + Parameters + ---------- + w : :class:`np.ndarray` [# obs wavelength bins] + Diagonal of matrix W. + M : :class:`np.ndarray` [# obs wavelength bins x # diags] + Resolution matrix, represented in sparse ROW form + giving nonzero entries in each row. + Jsp : :class:`tuple` + column-sparse matrix (endpts, J), where endpts = (s,e) + gives a half-open range of indices [s, e) with values + for this column, and P[:e-s] contains these values. + + Returns + ------- + Product WMJ in column-sparse form as tuple (endpts, P) + where endpts = (s,e) gives a half-open range of indices + [s, e) with values for this column, and P[:e-s] contains these + values. + + Note + ---- + We assume that the column-sparse input Jsp has enough space + allocated that we can overwrite each column of Jsp with the + corresponding column of P. + + """ + endpts, J = Jsp nbins, ndiag = M.shape ncol, maxColSize = J.shape @@ -598,22 +716,33 @@ def mulWMJ(w, M, jac): return (endpts, J) -# -# _prepare_bins -# Convert N bin centers to the info needed by the optimizer, -# returned as an opaque tuple. Given N bins, return values -# include -# -# - N+1 log bin edges. Edges are placed halfway between centers, -# with extrapolation at the ends. We return the natural log of -# each edge's wavelength, since that is what model and Jacobian -# building need. -# - N inverse bin widths, in an array padded by one cell on -# the left and right. -# -@jit(nopython=True, fastmath=False, nogil=True) + +@jit(nopython=True, nogil=True) def _prepare_bins(centers, camerapix): - + """ + Convert bin centers to the info needed by the optimizer, + returned as an opaque tuple. + + Parameters + ---------- + centers : :class:`np.ndarray` [# obs wavelength bins] + Center wavelength of each obs bin + camerapix : :class:`np.ndarray` of `int` [# cameras x 2] + pixels corresponding to each camera in obs wavelength array + + Returns + ------- + Tuple containing + - array of log bin edges. Edges are placed halfway between centers, + with extrapolation at the ends. We return the natural log of + each edge's wavelength, since that is what model and Jacobian + building need. + - array of inverse widths for each wavelength bin. The array is + zero-padded by one cell on the left and right to accomodate + edge-to-bin computations. + + """ + ncameras = camerapix.shape[0] edges = np.empty(len(centers) + ncameras, dtype=centers.dtype) ibin_widths = np.empty(len(centers) + 2, dtype=centers.dtype) From e150b0ca3bb1acaef87a9dda350c75587581402d Mon Sep 17 00:00:00 2001 From: jdbuhler Date: Tue, 3 Sep 2024 16:25:44 -0400 Subject: [PATCH 12/20] * update util.py for docstring convention * consistently import jit from numba instead of saying '@numba.jit' --- py/fastspecfit/continuum.py | 6 +-- py/fastspecfit/util.py | 86 +++++++++++++++++++++---------------- 2 files changed, 52 insertions(+), 40 deletions(-) diff --git a/py/fastspecfit/continuum.py b/py/fastspecfit/continuum.py index 8d18db9e..5b8e6fa9 100644 --- a/py/fastspecfit/continuum.py +++ b/py/fastspecfit/continuum.py @@ -7,7 +7,7 @@ """ import time import numpy as np -import numba +from numba import jit from fastspecfit.logger import log from fastspecfit.photometry import Photometry @@ -589,7 +589,7 @@ def call_nnls(modelflux, flux, ivar, xparam=None, debug=False, @staticmethod - @numba.jit(nopython=True, fastmath=True, nogil=True) + @jit(nopython=True, nogil=True, fastmath=True) def attenuate(M, A, zfactors, wave, dustflux): """ Compute attenuated version of a model spectrum, @@ -628,7 +628,7 @@ def attenuate(M, A, zfactors, wave, dustflux): M[i] = (M[i] + lbol_diff * dustflux[i]) * zfactors[i] @staticmethod - @numba.jit(nopython=True, fastmath=True, nogil=True) + @jit(nopython=True, nogil=True, fastmath=True) def attenuate_nodust(M, A, zfactors): """ Compute attenuated version of a model spectrum M, diff --git a/py/fastspecfit/util.py b/py/fastspecfit/util.py index e273d3c2..4aaba632 100644 --- a/py/fastspecfit/util.py +++ b/py/fastspecfit/util.py @@ -6,27 +6,29 @@ """ import numpy as np -import numba +from numba import jit from fastspecfit.logger import log -try: # this fails when building the documentation +try: # this fails when building the documentation from scipy import constants - C_LIGHT = constants.c / 1000.0 # [km/s] + C_LIGHT = constants.c / 1000.0 # [km/s] except: - C_LIGHT = 299792.458 # [km/s] + C_LIGHT = 299792.458 # [km/s] -FLUXNORM = 1e17 # flux normalization factor for all DESI spectra [erg/s/cm2/A] +FLUXNORM = 1e17 # flux normalization factor for all DESI spectra [erg/s/cm2/A] -# -# A BoxedScalar is an item of an Numpy -# structured scalar type that is initialized -# to all zeros and can then be passed -# around by reference. Access the .value -# field to unbox the scalar. -# class BoxedScalar(object): + """ + A BoxedScalar is an item of an Numpy + structured scalar type that is initialized + to all zeros and can then be passed + around by reference. Access the .value + field to unbox the scalar. + + """ + def __init__(self, dtype): self.value = np.zeros(1, dtype=dtype)[0] @@ -37,24 +39,28 @@ def __setitem__(self, key, v): self.value[key] = v -# -# A Pool encapsulates paaallel execution with a -# multiprocessing.Pool, falling back to sequential -# execution in the current process if just one worker -# is requested. -# -# Unlike multiprocessing.Pool, our starmap() function -# takes a list of keyword argument dictionaries, -# rather than a list of positional arguments. -# class MPPool(object): + """ + A Pool encapsulates paraallel execution with a + multiprocessing.Pool, falling back to sequential + execution in the current process if just one worker + is requested. - # create a pool with nworkers workers, using the current - # process if nworkers is 1. If initiializer is not None, - # apply this function to the arguments in keyword dictionary - # init_argdict on startup in each each worker subprocess. - def __init__(self, nworkers, initializer=None, init_argdict=None): + Unlike multiprocessing.Pool, our starmap() function + takes a list of keyword argument dictionaries, + rather than a list of positional arguments. + """ + + def __init__(self, nworkers, initializer=None, init_argdict=None): + """ + create a pool with nworkers workers, using the current + process if nworkers is 1. If initiializer is not None, + apply this function to the arguments in keyword dictionary + init_argdict on startup in each each worker subprocess. + + """ + initfunc = None if initializer is None else self.apply_to_dict if nworkers > 1: @@ -69,10 +75,14 @@ def __init__(self, nworkers, initializer=None, init_argdict=None): else: self.pool = None - # apply function func to each of a list of inputs, represented - # as a list of keyword argument dictionaries. + def starmap(self, func, argdicts): + """ + apply function func to each of a list of inputs, represented + as a list of keyword argument dictionaries. + """ + # we cannot pickle a local function, so we must pass # both func and the argument dictionary to the subprocess # worker and have it apply one to the other. @@ -86,12 +96,14 @@ def starmap(self, func, argdicts): return out - # close our multiprocess pool if we created one def close(self): + """ + close our multiprocess pool if we created one + """ + if self.pool is not None: self.pool.close() - @staticmethod def apply_to_dict(f, argdict): return f(**argdict) @@ -335,7 +347,7 @@ def minfit(x, y, return_coeff=False): # array copies and redundant summation # on each iteration # -@numba.jit(nopython=True, nogil=True) +@jit(nopython=True, nogil=True) def sigmaclip(c, low=3., high=3.): n = len(c) @@ -372,19 +384,19 @@ def sigmaclip(c, low=3., high=3.): # Numba's quantile impl is much faster # than Numpy's standard version -@numba.jit(nopython=True, nogil=True) +@jit(nopython=True, nogil=True) def quantile(A, q): return np.quantile(A, q) # Numba's median impl is also faster -@numba.jit(nopython=True, nogil=True) +@jit(nopython=True, nogil=True) def median(A): return np.median(A) # Open-coded Numba trapz is much faster than np.traz -@numba.jit(nopython=True, fastmath=True, nogil=True) +@jit(nopython=True, nogil=True, fastmath=True) def trapz(y, x): res = 0. for i in range(len(x) - 1): @@ -453,7 +465,7 @@ def trapz_rebin_pre(bin_centers): return (edges, ibw) -@numba.jit(nopython=True, fastmath=True, nogil=True) +@jit(nopython=True, nogil=True, fastmath=True) def _trapz_rebin(x, y, edges, ibw, out): """ Trapezoidal rebinning @@ -519,7 +531,7 @@ def y_at(edge_x, j): # j: largest j s.t. x[j] < edge_x return results -@numba.jit(nopython=True, nogil=True) +@jit(nopython=True, nogil=True) def centers2edges(centers): """ Convert bin centers to bin edges, guessing at what you probably meant. From c1e5d2ea3cdb927e6036902111f21ad73ad7fb6d Mon Sep 17 00:00:00 2001 From: jdbuhler Date: Wed, 4 Sep 2024 11:23:24 -0400 Subject: [PATCH 13/20] remove stray character before comment --- py/fastspecfit/igm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/py/fastspecfit/igm.py b/py/fastspecfit/igm.py index 3dad887d..2efe9c96 100644 --- a/py/fastspecfit/igm.py +++ b/py/fastspecfit/igm.py @@ -11,7 +11,7 @@ class Inoue14(object): - r""" + """ IGM absorption from Inoue et al. (2014) Parameters From 23aae6f3c9e420e2d5d889fa0344af675eba0677 Mon Sep 17 00:00:00 2001 From: jdbuhler Date: Wed, 4 Sep 2024 11:29:58 -0400 Subject: [PATCH 14/20] fix docstrings to conform to rest of code --- py/fastspecfit/emlines.py | 76 ++++++++++++++++----------------------- 1 file changed, 30 insertions(+), 46 deletions(-) diff --git a/py/fastspecfit/emlines.py b/py/fastspecfit/emlines.py index 3ca13057..90e7ccfa 100644 --- a/py/fastspecfit/emlines.py +++ b/py/fastspecfit/emlines.py @@ -124,8 +124,9 @@ def __init__(self, emline_table, uniqueid=None, stronglines=False): self.param_table['modelname'] = \ np.array([ s.replace('_amp', '_modelamp').upper() for s in param_names ]) - # Record which lines are within the limits of the cameras + def compute_inrange_lines(self, redshift, wavelims=(3600., 9900.), wavepad=5*0.8): + """ Record which lines are within the limits of the cameras """ zlinewave = self.line_table['restwave'].value * (1. + redshift) @@ -133,27 +134,33 @@ def compute_inrange_lines(self, redshift, wavelims=(3600., 9900.), wavepad=5*0.8 ((zlinewave > (wavelims[0] + wavepad)) & \ (zlinewave < (wavelims[1] - wavepad))) - # Build emission line model tables, with and without suppression of broad lines. - # Establish fixed (i.e., forced to zero) params, tying relationships between - # params, and doublet relationships for each model, as well as the relationship - # between lines and their parameters, which we record in the line table. - # - # Parameter fixing needs to know which lines are within the observed wavelength - # ranges of the cameras, so we first add this information to the line table. + def build_linemodels(self, separate_oiii_fit=True): - # create_model() - # - # Given the tying info for the model and the list of in-range - # lines, determine which parameters of the model are fixed and - # free, and create the model's table. - # - # We fix all parameters for out-of-range lines to zero, unless - # the parameter is tied to a parameter for an in-range line. - # We also allow the caller to specify parameters that should - # be fixed regardless. + """Build emission line model tables, with and without + suppression of broad lines. Establish fixed (i.e., forced to + zero) params, tying relationships between params, and doublet + relationships for each model, as well as the relationship + between lines and their parameters, which we record in the + line table. + + Parameter fixing needs to know which lines are within the + observed wavelength ranges of the cameras, so we first add + this information to the line table. + """ + def create_model(tying_info, forceFixed=[]): + """Given the tying info for the model and the list of + in-range lines, determine which parameters of the model + are fixed and free, and create the model's table. + + We fix all parameters for out-of-range lines to zero, + unless the parameter is tied to a parameter for an + in-range line. We also allow the caller to specify + parameters that should be fixed regardless. + """ + n_params = len(self.param_table) isfixed = np.full(n_params, False, dtype=bool) @@ -346,11 +353,9 @@ def tie_line(tying_info, line_params, source_linename, amp_factor=None): return linemodel_broad, linemodel_nobroad - # debugging function to print basic contents of a line model def summarize_linemodel(self, linemodel): - """Simple function to summarize an input linemodel.""" - + def _print(line_mask): for line in np.where(line_mask)[0]: line_name = self.line_table['name'][line] @@ -601,27 +606,6 @@ def optimize(self, linemodel, # least_squares wants two arrays, not a 2D array bounds = ( bounds[:,0], bounds[:,1] ) - - """ - ######################## - # SANITY CHECK - ######################## - if continuum_patches is not None: - J = jac(initial_guesses) - delta = 1e-5 - - v = np.zeros(len(initial_guesses)) - for p in range(len(v)): - v[p] = 1. - ddp = (objective(initial_guesses + delta*v) - objective(initial_guesses - delta*v))/(2*delta) - ddpj = J.dot(v) - v[p] = 0. - - d = np.abs(ddpj - ddp) - if d > 1e-5: - print("JAC ", p, d) - """ - fit_info = least_squares(objective, initial_guesses, jac=jac, args=(), max_nfev=5000, xtol=1e-10, ftol=1e-5, #x_scale='jac' gtol=1e-10, tr_solver='lsmr', tr_options={'maxiter': 1000, 'regularize': True}, @@ -1301,12 +1285,12 @@ def emline_specfit(data, result, continuummodel, smooth_continuum, result['RCHI2_LINE'] = finalchi2 #result['NDOF_LINE'] = finalndof - result['DELTA_LINECHI2'] = delta_linechi2_balmer # chi2_nobroad - chi2_broad - result['DELTA_LINENDOF'] = delta_linendof_balmer # ndof_nobroad - ndof_broad + result['DELTA_LINECHI2'] = delta_linechi2_balmer # chi2_nobroad - chi2_broad + result['DELTA_LINENDOF'] = delta_linendof_balmer # ndof_nobroad - ndof_broad # full-fit reduced chi2 rchi2 = np.sum(oemlineivar * (specflux - (continuummodelflux + smoothcontinuummodelflux + emmodel))**2) - rchi2 /= np.sum(oemlineivar > 0) # dof?? + rchi2 /= np.sum(oemlineivar > 0) # dof?? result['RCHI2'] = rchi2 # I believe that all the elements of the coadd_wave vector are contained @@ -1435,7 +1419,7 @@ def synthphot_spectrum(phot, data, result, modelwave, modelflux): filters = phot.synth_filters[data['photsys']] synthmaggies = Photometry.get_ab_maggies(filters, modelflux / FLUXNORM, modelwave) - model_synthmag = Photometry.to_nanomaggies(synthmaggies) # units of nanomaggies + model_synthmag = Photometry.to_nanomaggies(synthmaggies) # units of nanomaggies model_synthphot = Photometry.parse_photometry(phot.synth_bands, maggies=synthmaggies, nanomaggies=False, From 87d4edbb0def42992ae53e6f9a8702bcd2d187ff Mon Sep 17 00:00:00 2001 From: jdbuhler Date: Thu, 5 Sep 2024 08:50:50 -0400 Subject: [PATCH 15/20] * revert removal of raw tag from string in igm.py * make log msg about using mkl_fft debug instead of info * minor whitespace cleanups --- py/fastspecfit/continuum.py | 4 ++-- py/fastspecfit/igm.py | 2 +- py/fastspecfit/templates.py | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/py/fastspecfit/continuum.py b/py/fastspecfit/continuum.py index 5b8e6fa9..c3525c41 100644 --- a/py/fastspecfit/continuum.py +++ b/py/fastspecfit/continuum.py @@ -751,9 +751,9 @@ def continuum_to_spectroscopy(self, contmodel): camerapix = self.data['camerapix'] specwave = self.data['wave'] specres = self.data['res'] - + modelflux = np.empty(self.wavelen) - + for icam, (s, e) in enumerate(camerapix): resampflux = trapz_rebin(self.ztemplatewave, contmodel, diff --git a/py/fastspecfit/igm.py b/py/fastspecfit/igm.py index 2efe9c96..3dad887d 100644 --- a/py/fastspecfit/igm.py +++ b/py/fastspecfit/igm.py @@ -11,7 +11,7 @@ class Inoue14(object): - """ + r""" IGM absorption from Inoue et al. (2014) Parameters diff --git a/py/fastspecfit/templates.py b/py/fastspecfit/templates.py index cf461de8..72d3d9af 100644 --- a/py/fastspecfit/templates.py +++ b/py/fastspecfit/templates.py @@ -206,9 +206,9 @@ def init_ffts(self): if find_spec("mkl_fft") is not None: import mkl_fft._scipy_fft_backend as be sc_fft.set_global_backend(be) - + self.convolve = sc_sig.convolve - log.info('Using mkl_fft library for FFTs') + log.debug('Using mkl_fft library for FFTs') else: self.convolve = sc_sig.oaconvolve From 9db05310f7b487ff040de84209fbb0c7c2d9acb4 Mon Sep 17 00:00:00 2001 From: jdbuhler Date: Thu, 5 Sep 2024 13:59:28 -0400 Subject: [PATCH 16/20] remove trailing whitespace from all lines, per house style --- py/fastspecfit/emline_fit/interface.py | 160 ++++++++++---------- py/fastspecfit/emline_fit/jacobian.py | 74 ++++----- py/fastspecfit/emline_fit/model.py | 80 +++++----- py/fastspecfit/emline_fit/params_mapping.py | 88 +++++------ py/fastspecfit/emline_fit/sparse_rep.py | 76 +++++----- py/fastspecfit/emline_fit/utils.py | 14 +- 6 files changed, 246 insertions(+), 246 deletions(-) diff --git a/py/fastspecfit/emline_fit/interface.py b/py/fastspecfit/emline_fit/interface.py index 8759d4b1..4ff86c4a 100644 --- a/py/fastspecfit/emline_fit/interface.py +++ b/py/fastspecfit/emline_fit/interface.py @@ -28,7 +28,7 @@ class EMLine_Objective(object): """ Objective function for emission-line fitting. """ - + def __init__(self, obs_bin_centers, obs_fluxes, @@ -61,9 +61,9 @@ def __init__(self, else None. """ - + self.dtype = obs_fluxes.dtype - + self.obs_fluxes = obs_fluxes self.obs_weights = obs_weights self.redshift = redshift @@ -78,7 +78,7 @@ def __init__(self, self.patch_endpts = continuum_patches['endpts'].value self.patch_pivotwave = continuum_patches['pivotwave'].value - + self.J_P = patch_jacobian(obs_bin_centers, obs_weights, self.patch_endpts, @@ -86,10 +86,10 @@ def __init__(self, else: self.nPatches = 0 self.J_P = None - + self.log_obs_bin_edges, self.ibin_widths = \ _prepare_bins(obs_bin_centers, camerapix) - + # temporary storage to prevent allocation in params_mapping # on every call to objective/jacobian self.line_parameters = \ @@ -111,12 +111,12 @@ def objective(self, free_parameters): Returns ------- `np.ndarray` of residuals for each wavelength bin. - + """ - + nLineParameters = len(free_parameters) - 2 * self.nPatches line_free_parameters = free_parameters[:nLineParameters] - + # # expand line free parameters into complete # line parameter array, handling tied params @@ -124,9 +124,9 @@ def objective(self, free_parameters): # line_parameters = self.params_mapping.mapFreeToFull(line_free_parameters, out=self.line_parameters) - + model_fluxes = np.empty_like(self.obs_fluxes, dtype=self.dtype) - + _build_model_core(line_parameters, self.line_wavelengths, self.redshift, @@ -146,13 +146,13 @@ def objective(self, free_parameters): self.patch_pivotwave, slopes, intercepts) - + # turn model fluxes into residuals in-place to avoid # unwanted memory allocation residuals = model_fluxes residuals -= self.obs_fluxes residuals *= self.obs_weights - + return residuals @@ -172,18 +172,18 @@ def jacobian(self, free_parameters): ------- Sparse Jacobian at given parameter value as defined in sparse_rep.py. - + """ - + nLineFreeParms = len(free_parameters) - 2 * self.nPatches line_free_parameters = free_parameters[:nLineFreeParms] - + # # expand free paramters into complete # parameter array, handling tied params # and doublets # - + line_parameters = self.params_mapping.mapFreeToFull(line_free_parameters, out=self.line_parameters) @@ -198,7 +198,7 @@ def jacobian(self, free_parameters): # exists to either side of this range, so we can pass # those in as dummies. ibw = self.ibin_widths[s:e+3] - + idealJac = \ emline_model_jacobian(line_parameters, self.log_obs_bin_edges[s+icam:e+icam+1], @@ -206,15 +206,15 @@ def jacobian(self, free_parameters): self.redshift, self.line_wavelengths, self.resolution_matrices[icam].ndiag) - + # ignore any columns corresponding to fixed parameters endpts = idealJac[0] endpts[self.params_mapping.fixedMask(), :] = 0 - + jacs.append( mulWMJ(self.obs_weights[s:e], self.resolution_matrices[icam].rowdata(), idealJac) ) - + nBins = np.sum(np.diff(self.camerapix)) nFreeParms = len(free_parameters) J = EMLineJacobian((nBins, nFreeParms), nLineFreeParms, @@ -236,7 +236,7 @@ def build_model(redshift, continuum_patches=None): """ Compatibility entry point to compute modeled fluxes. - + Parameters ---------- redshift : :class:`np.float64` @@ -256,11 +256,11 @@ def build_model(redshift, else None. """ - + log_obs_bin_edges, ibin_widths = _prepare_bins(obs_bin_centers, camerapix) - + model_fluxes = np.empty_like(obs_bin_centers, dtype=obs_bin_centers.dtype) - + _build_model_core(line_parameters, line_wavelengths, redshift, @@ -272,15 +272,15 @@ def build_model(redshift, # suppress negative pixels arising from resolution matrix model_fluxes[model_fluxes < 0.] = 0. - + if continuum_patches is not None: # add patch pedestals to model fluxes _add_patches(obs_bin_centers, model_fluxes, continuum_patches['endpts'].value, - continuum_patches['pivotwave'].value, + continuum_patches['pivotwave'].value, continuum_patches['slope'].value, continuum_patches['intercept'].value) - + return model_fluxes @@ -291,7 +291,7 @@ class MultiLines(object): more cameras. Return an object that lets caller obtain a rendered model for each individual line as a sparse array. """ - + def __init__(self, line_parameters, obs_bin_centers, @@ -319,7 +319,7 @@ def __init__(self, Pixels corresponding to each camera in obs wavelength array. """ - + @jit(nopython=True, nogil=True) def _suppress_negative_fluxes(endpts, M): """ @@ -329,7 +329,7 @@ def _suppress_negative_fluxes(endpts, M): s, e = endpts[i] for j in range(e-s): M[i,j] = np.maximum(M[i,j], 0.) - + self.line_models = [] _build_multimodel_core(line_parameters, obs_bin_centers, @@ -364,7 +364,7 @@ def getLine(self, line): the bin values. """ - + s = 1000000000 e =-1 @@ -374,14 +374,14 @@ def getLine(self, line): # observed flux array. live_models = [] for i, line_model in enumerate(self.line_models): - + ls, le = line_model[0][line] - + if ls < le: # line has nonzero flux bins s = np.minimum(s, ls) e = np.maximum(e, le) live_models.append(i) - + if len(live_models) == 0: # line has no nonzero flux bins return (0, 0), np.empty((0), dtype=np.float64) @@ -402,7 +402,7 @@ def getLine(self, line): ldata = self.line_models[i][1][line] data[ls-s:le-s] = ldata[:le-ls] - + return (s, e), data @@ -437,7 +437,7 @@ def find_peak_amplitudes(line_parameters, Array with maximum amplitude observed for each line. """ - + @jit(nopython=True, nogil=True) def _update_line_maxima(max_amps, line_models): """ @@ -447,7 +447,7 @@ def _update_line_maxima(max_amps, line_models): in place. """ endpts, vals = line_models - + # find the highest flux for each peak; if it's # bigger than any seen so far, update global max for i in range(vals.shape[0]): @@ -465,10 +465,10 @@ def _update_line_maxima(max_amps, line_models): resolution_matrices, camerapix, lambda m: _update_line_maxima(max_amps, m)) - + return max_amps - + ########################################################################## @@ -483,7 +483,7 @@ def _build_model_core(line_parameters, """ Core loop for computing a combined model flux from a set of spectral emission lines. - + Parameters ---------- line_parameters : :class:`np.ndarray` @@ -504,9 +504,9 @@ def _build_model_core(line_parameters, Returns computed model flux for each wavelength bin. """ - + for icam, campix in enumerate(camerapix): - + # start and end for obs fluxes of camera icam s, e = campix @@ -515,13 +515,13 @@ def _build_model_core(line_parameters, # exists to either side of this range, so we can pass # those in as dummies. ibw = ibin_widths[s:e+3] - + mf = emline_model(line_wavelengths, line_parameters, log_obs_bin_edges[s+icam:e+icam+1], redshift, ibw) - + # convolve model with resolution matrix and store in # this camera's subrange of model_fluxes resolution_matrices[icam].dot(mf, model_fluxes[s:e]) @@ -538,7 +538,7 @@ def _build_multimodel_core(line_parameters, Core loop for computing array of individual line flux models sparsely from a set of spectral emission lines. Result is not returned but passed to a consumer function supplied by the caller. - + Parameters ---------- line_parameters : :class:`np.ndarray` @@ -558,21 +558,21 @@ def _build_multimodel_core(line_parameters, for each camera. """ - + log_obs_bin_edges, ibin_widths = _prepare_bins(obs_bin_centers, camerapix) - + for icam, campix in enumerate(camerapix): - + # start and end for obs fluxes of camera icam s, e = campix - + # Actual inverse bin widths are in ibin_widths[s+1:e+2]. # Setup guarantees that at least one more array entry # exists to either side of this range, so we can pass # those in as dummies. ibw = ibin_widths[s:e+3] - + # compute model waveform for each spectral line line_models = emline_perline_models(line_wavelengths, line_parameters, @@ -580,15 +580,15 @@ def _build_multimodel_core(line_parameters, redshift, ibw, resolution_matrices[icam].ndiag) - + # convolve each line's waveform with resolution matrix endpts, M = mulWMJ(np.ones(e - s), resolution_matrices[icam].rowdata(), line_models) - + # adjust endpoints to reflect camera range endpts += s - + consumer_fun((endpts, M)) @@ -620,21 +620,21 @@ def _add_patches(obs_bin_centers, Intercept of each patch pedestal. """ - + # add patch pedestals to line model nPatches = len(slopes) - + for ipatch in range(nPatches): s, e = patch_endpts[ipatch] pivotwave = patch_pivotwaves[ipatch] - + slope = slopes[ipatch] intercept = intercepts[ipatch] - + for j in range(s,e): model_fluxes[j] += \ slope * (obs_bin_centers[j] - pivotwave) + intercept - + ################################################################### @@ -645,7 +645,7 @@ def mulWMJ(w, M, Jsp): W is a diagonal weight matrix M is a resolution matrix J is a column-sparse matrix giving the nonzero - values in one contiguous range per column + values in one contiguous range per column Parameters ---------- @@ -656,7 +656,7 @@ def mulWMJ(w, M, Jsp): giving nonzero entries in each row. Jsp : :class:`tuple` column-sparse matrix (endpts, J), where endpts = (s,e) - gives a half-open range of indices [s, e) with values + gives a half-open range of indices [s, e) with values for this column, and P[:e-s] contains these values. Returns @@ -674,35 +674,35 @@ def mulWMJ(w, M, Jsp): """ endpts, J = Jsp - + nbins, ndiag = M.shape ncol, maxColSize = J.shape - + hdiag = ndiag//2 - + # temporary buffer for each column of WMJ buf = np.empty(maxColSize, dtype=J.dtype) - + for j in range(ncol): # boundaries of nonzero entries # in jth column of J s, e = endpts[j] - + if s == e: # no nonzero values in column j continue - + # boundaries of entries in jth column of P # impacted by matrix multiply imin = np.maximum(s - hdiag, 0) imax = np.minimum(e + hdiag, nbins) # one past last impacted entry - + for i in range(imin, imax): - + # boundaries of interval of k where both # M[i, k] and J[k, j] are nonzero. kmin = np.maximum(i - hdiag, s) kmax = np.minimum(i + hdiag, e - 1) - + acc = 0. for k in range(kmin, kmax + 1): acc += M[i, k - i + hdiag] * J[j, k - s] @@ -713,7 +713,7 @@ def mulWMJ(w, M, Jsp): newE = np.minimum(imax, nbins) J[j, :newE - newS] = buf[:newE - newS] endpts[j] = np.array([newS, newE]) - + return (endpts, J) @@ -739,36 +739,36 @@ def _prepare_bins(centers, camerapix): building need. - array of inverse widths for each wavelength bin. The array is zero-padded by one cell on the left and right to accomodate - edge-to-bin computations. + edge-to-bin computations. """ - + ncameras = camerapix.shape[0] edges = np.empty(len(centers) + ncameras, dtype=centers.dtype) ibin_widths = np.empty(len(centers) + 2, dtype=centers.dtype) - + for icam, campix in enumerate(camerapix): - + s, e = campix icenters = centers[s:e] - + #- interior edges are just points half way between bin centers int_edges = 0.5 * (icenters[:-1] + icenters[1:]) - + #- exterior edges are extrapolation of interior bin sizes edge_l = icenters[ 0] - (icenters[ 1] - int_edges[ 0]) edge_r = icenters[-1] + (icenters[-1] - int_edges[-1]) - + edges[s + icam] = edge_l edges[s + icam + 1:e + icam] = int_edges edges[e + icam] = edge_r # add 1 to indices i ibin_widths to skip dummy at 0 ibin_widths[s+1:e+1] = 1. / np.diff(edges[s+icam : e+icam+1]) - + # dummies before and after widths are needed # for corner cases in edge -> bin computation ibin_widths[0] = 0. ibin_widths[-1] = 0. - + return (np.log(edges), ibin_widths) diff --git a/py/fastspecfit/emline_fit/jacobian.py b/py/fastspecfit/emline_fit/jacobian.py index 55e289a9..741fdb13 100644 --- a/py/fastspecfit/emline_fit/jacobian.py +++ b/py/fastspecfit/emline_fit/jacobian.py @@ -47,55 +47,55 @@ def emline_model_jacobian(line_parameters, """ SQRT_2PI = np.sqrt(2*np.pi) - + nbins = len(log_obs_bin_edges) - 1 - + line_amplitudes, line_vshifts, line_sigmas = \ np.split(line_parameters, 3) - + # buffers for per-parameter calculations, sized large # enough for max possible range [s .. e), plus extra # padding as directed by caller. max_width = max_buffer_width(log_obs_bin_edges, line_sigmas, padding) - + nlines = len(line_wavelengths) dd = np.empty((3 * nlines, max_width), dtype=line_amplitudes.dtype) endpts = np.zeros((3 * nlines, 2), dtype=np.int32) starts = endpts[:, 0] ends = endpts[:, 1] - + # compute partial derivatives for avg values of all Gaussians # inside each bin. For each Gaussian, we only compute # contributions for bins where it is non-negligible. for j in range(len(line_wavelengths)): - + # line width sigma = line_sigmas[j] / C_LIGHT - + c0 = SQRT_2PI * np.exp(0.5 * sigma**2) - + # wavelength shift for spectral lines line_shift = 1. + redshift + line_vshifts[j] / C_LIGHT shifted_line = line_wavelengths[j] * line_shift log_shifted_line = np.log(shifted_line) - - # leftmost edge i that needs a value (> 0) for this line + + # leftmost edge i that needs a value (> 0) for this line lo = np.searchsorted(log_obs_bin_edges, log_shifted_line - MAX_SDEV * sigma, side="left") - + # leftmost edge i that does *not* need a value (== 1) for this line hi = np.searchsorted(log_obs_bin_edges, log_shifted_line + MAX_SDEV * sigma, side="right") # check if entire Gaussian is outside bounds of log_obs_bin_edges - if hi == 0 or lo == len(log_obs_bin_edges): + if hi == 0 or lo == len(log_obs_bin_edges): continue - + nedges = hi - lo + 2 # compute values at edges [lo - 1 ... hi] - + # Compute contribs of each line to each partial derivative in place. # No sharing of params between peaks means that we never have to # add contributions from two peaks to same line. @@ -104,34 +104,34 @@ def emline_model_jacobian(line_parameters, dds_vals = dd[2*nlines + j] offset = (log_shifted_line / sigma + sigma) if sigma > 0. else 0. - + c = c0 * line_wavelengths[j] A = c / C_LIGHT * line_amplitudes[j] - + # vals[i] --> edge i + lo - 1 - + dda_vals[0] = 0. # edge lo - 1 ddv_vals[0] = 0. dds_vals[0] = 0. - + for i in range(1, nedges - 1): - + # x - offset = (log(lambda_j) - mu_i)/sigma - sigma, # the argument of the Gaussian integral - + x = log_obs_bin_edges[i+lo-1]/sigma - offset pdf = norm_pdf(x) cdf = norm_cdf(x) - + dda_vals[i] = c * line_shift * sigma * cdf ddv_vals[i] = A * (sigma * cdf - pdf) dds_vals[i] = A * line_shift * \ ((1 + sigma**2) * cdf - (x + 2*sigma) * pdf) - + dda_vals[nedges - 1] = c * line_shift * sigma # edge hi ddv_vals[nedges - 1] = A * sigma dds_vals[nedges - 1] = A * line_shift * (1 + sigma**2) - + # Compute partial derivs for bin i+lo-1 for 0 <= i < nedges - 1. # But: # * if lo == 0, bin lo - 1 is not defined, so skip it @@ -145,7 +145,7 @@ def emline_model_jacobian(line_parameters, dlo = 1 if lo == 0 else 0 dhi = 1 if hi == len(log_obs_bin_edges) else 0 - + for i in range(dlo, nedges - 1 - dhi): dda_vals[i-dlo] = (dda_vals[i+1] - dda_vals[i]) * ibin_widths[i+lo] ddv_vals[i-dlo] = (ddv_vals[i+1] - ddv_vals[i]) * ibin_widths[i+lo] @@ -153,17 +153,17 @@ def emline_model_jacobian(line_parameters, # starts[j] is set to first valid bin starts[j] = lo - 1 + dlo - + # ends[j] is set one past last valid bin ends[j] = hi - dhi - + # replicate first third of endpts (which is what we # set above) twice more, since same endpts apply to # all three params of each line for i in range(1,3): endpts[i*nlines:(i+1)*nlines,:] = endpts[:nlines,:] - # for lines with zero amplitude, + # for lines with zero amplitude, # partial derivatives w/r to their vshifts # and sigmas are zero for i, amp in enumerate(line_amplitudes): @@ -176,8 +176,8 @@ def emline_model_jacobian(line_parameters, for i, sig in enumerate(line_sigmas): if sig == 0.: endpts[i, :] = 0 - - return (endpts, dd) + + return (endpts, dd) @staticmethod @@ -203,7 +203,7 @@ def patch_jacobian(obs_bin_centers, Endpoints of each patch in wavelength array. patch_pivotwave : :class:`np.ndarray` [# patches] Wavelength offset for fitted affine params of each patch . - + Returns ------- :class:`tuple` (endpts, M) @@ -211,7 +211,7 @@ def patch_jacobian(obs_bin_centers, [ endpts[j,0] , endpts[j,1] ], which are stored in M[j]. """ - + nPatches = patch_endpts.shape[0] # @@ -219,16 +219,16 @@ def patch_jacobian(obs_bin_centers, # for slopes and intercepts) and compute the # maximum width of any patch. # - + endpts = np.empty((2*nPatches, 2), dtype=np.int32) maxPatchWidth = 0 for i in range(nPatches): s, e = patch_endpts[i] - + endpts[i] = (s, e) endpts[i+nPatches] = (s, e) - + maxPatchWidth = np.maximum(maxPatchWidth, e - s) # @@ -237,7 +237,7 @@ def patch_jacobian(obs_bin_centers, # patch. These derivatives are nonzero only # within the boundaries of the patch. # - + M = np.empty((2*nPatches, maxPatchWidth)) for i in range(nPatches): s, e = endpts[i] @@ -246,8 +246,8 @@ def patch_jacobian(obs_bin_centers, M[i,:e-s] = \ (obs_bin_centers[s:e] - patch_pivotwave[i]) * \ obs_weights[s:e] - + # dobj/dintercept for patch M[i + nPatches, :e-s] = obs_weights[s:e] # 1. x obs_weights - + return (endpts, M) diff --git a/py/fastspecfit/emline_fit/model.py b/py/fastspecfit/emline_fit/model.py index baf2006a..46d35845 100644 --- a/py/fastspecfit/emline_fit/model.py +++ b/py/fastspecfit/emline_fit/model.py @@ -38,14 +38,14 @@ def emline_model(line_wavelengths, Returns ------- - `np.ndarray` of length [# obs wavelength bins] with + `np.ndarray` of length [# obs wavelength bins] with average fluxes in each observed wavelength bin """ line_amplitudes, line_vshifts, line_sigmas = \ np.split(line_parameters, 3) - + # output per-bin fluxes nbins = len(log_obs_bin_edges) - 1 model_fluxes = np.zeros(nbins, dtype=line_amplitudes.dtype) @@ -58,10 +58,10 @@ def emline_model(line_wavelengths, # For each Gaussian, we only compute area contributions # for bins where it is non-negligible. for j in range(len(line_wavelengths)): - + if line_amplitudes[j] == 0.: continue - + s, e = emline_model_core(line_wavelengths[j], line_amplitudes[j], line_vshifts[j], @@ -70,12 +70,12 @@ def emline_model(line_wavelengths, redshift, ibin_widths, bin_vals) - + # add bin avgs for this peak to bins [s, e) model_fluxes[s:e] += bin_vals[:e-s] - + return model_fluxes - + @jit(nopython=True, nogil=True) def emline_perline_models(line_wavelengths, @@ -89,7 +89,7 @@ def emline_perline_models(line_wavelengths, line individually its waveform in the range of bins described by log_bin_edges. The waveforms are returned sparsely in the same forma as the Jacobian below. - + This function shares the core computation of emline_model() but does not collapse the lines into one composite. @@ -107,7 +107,7 @@ def emline_perline_models(line_wavelengths, Inverse widths of each observed wavelength bin. padding : :class:`int` Number of entries by which to pad each sparse row for later use. - + Returns ------- Row-sparse matrix of size [# lines x # obs wavelength bins] @@ -117,21 +117,21 @@ def emline_perline_models(line_wavelengths, (and may have some empty cells thereafter if padding is non-zero). """ - + line_amplitudes, line_vshifts, line_sigmas = \ np.split(line_parameters, 3) - + nbins = len(log_obs_bin_edges) - 1 - + # buffers for per-line calculations, sized large # enough for max possible range [s .. e), plus extra # padding as directed by caller. max_width = max_buffer_width(log_obs_bin_edges, line_sigmas, padding) - + nlines = len(line_wavelengths) line_profiles = np.empty((nlines, max_width), dtype=line_amplitudes.dtype) endpts = np.zeros((nlines, 2), dtype=np.int32) - + # compute total area of all Gaussians inside each bin. # For each Gaussian, we only compute area contributions # for bins where it is non-negligible. @@ -139,9 +139,9 @@ def emline_perline_models(line_wavelengths, if line_amplitudes[j] == 0.: continue - + bin_vals = line_profiles[j] - + s, e = emline_model_core(line_wavelengths[j], line_amplitudes[j], line_vshifts[j], @@ -153,7 +153,7 @@ def emline_perline_models(line_wavelengths, endpts[j,0] = s endpts[j,1] = e - + return (endpts, line_profiles) @@ -169,7 +169,7 @@ def emline_model_core(line_wavelength, """ Compute the flux contribution of one spectral line to a model spectrum. We compute the range of bins where the line - contributes flux, then write those contributions to the + contributes flux, then write those contributions to the output array vals and return half-open range of bins [s, e) to which it contributes. @@ -201,57 +201,57 @@ def emline_model_core(line_wavelength, Notes ----- vals MUST have length at least two more than the largest number - of bins that could be computed. Entries in vals after the returned + of bins that could be computed. Entries in vals after the returned range may be set to arbitrary values. """ - + SQRT_2PI = np.sqrt(2 * np.pi) - + # line width sigma = line_sigma / C_LIGHT - + c = SQRT_2PI * sigma * np.exp(0.5 * sigma**2) - + # wavelength shift for spectral lines line_shift = 1. + redshift + line_vshift / C_LIGHT - + shifted_line = line_wavelength * line_shift log_shifted_line = np.log(shifted_line) - - # leftmost edge i that needs a value (> 0) for this line + + # leftmost edge i that needs a value (> 0) for this line lo = np.searchsorted(log_obs_bin_edges, log_shifted_line - MAX_SDEV * sigma, side="left") - + # leftmost edge i that does *not* need a value (== 1) for this line hi = np.searchsorted(log_obs_bin_edges, log_shifted_line + MAX_SDEV * sigma, side="right") # check if entire Gaussian is outside bounds of log_obs_bin_edges - if hi == 0 or lo == len(log_obs_bin_edges): + if hi == 0 or lo == len(log_obs_bin_edges): return (0,0) - + nedges = hi - lo + 2 # compute values at edges [lo - 1 ... hi] - + A = c * line_amplitude * shifted_line offset = (log_shifted_line / sigma + sigma) if sigma > 0. else 0. - + # vals[i] --> edge i + lo - 1 - + vals[0] = 0. # edge lo - 1 - + for i in range(1, nedges-1): - + # x = (log(lambda_j) - mu_i)/sigma - sigma, # the argument of the Gaussian integral - + x = log_obs_bin_edges[i+lo-1]/sigma - offset vals[i] = A * norm_cdf(x) - + vals[nedges-1] = A # edge hi - + # Compute avg of bin i+lo-1 for 0 <= i < nedges - 1. # But: # * if lo == 0, bin lo - 1 is not defined, so skip it @@ -266,11 +266,11 @@ def emline_model_core(line_wavelength, # We ensure that values for all valid bins are written # starting at vals[0], and return the range of indices # of the computed bins w/r to the input bin array. - + dlo = 1 if lo == 0 else 0 dhi = 1 if hi == len(log_obs_bin_edges) else 0 - + for i in range(dlo, nedges - 1 - dhi): vals[i-dlo] = (vals[i+1] - vals[i]) * ibin_widths[i+lo] - + return (lo - 1 + dlo, hi - dhi) diff --git a/py/fastspecfit/emline_fit/params_mapping.py b/py/fastspecfit/emline_fit/params_mapping.py index 527ef15c..c04ac9d2 100644 --- a/py/fastspecfit/emline_fit/params_mapping.py +++ b/py/fastspecfit/emline_fit/params_mapping.py @@ -11,7 +11,7 @@ class ParamsMapping(object): work to update it for each new set of free parameters. """ - + def __init__(self, nParms, isFree, tiedSources, tiedFactors, @@ -33,12 +33,12 @@ def __init__(self, nParms, Parameters that are doublet ratios vs. a source. doubletSources: :class:`np.ndarray` of `int` Source parameters for all doublet ratios. - + """ - + self.nParms = nParms self.nFreeParms = np.sum(isFree) - + # permutation mapping each free parameter in full list # to its location in free list pFree = np.empty(self.nParms, dtype=np.int32) @@ -47,17 +47,17 @@ def __init__(self, nParms, self._precomputeMapping(isFree, pFree, tiedSources, tiedFactors, doubletTargets, doubletSources) - + self._precomputeJacobian() - + def fixedMask(self): """ Return Boolean mask with "True" for each fixed parameter. """ return self.isFixed - + def mapFreeToFull(self, freeParms, out=None, patchDoublets=True): """ Given a vector of free parameters, return the corresponding @@ -90,7 +90,7 @@ def mapFreeToFull(self, freeParms, out=None, patchDoublets=True): self.doubletPatches, out, patchDoublets) - + @staticmethod @jit(nopython=True, nogil=True) @@ -102,15 +102,15 @@ def _mapFreeToFull(freeParms, nParms, sources, factors, if fullParms is None: fullParms = np.empty(nParms, dtype=freeParms.dtype) - + for j, src_j_free in enumerate(sources): fullParms[j] = factors[j] # copy fixed zeros if src_j_free != -1: fullParms[j] *= freeParms[src_j_free] - + return fullParms - - + + def getJacobian(self, freeParms): """ Given a vector v of free parameters, return the Jacobian @@ -127,14 +127,14 @@ def getJacobian(self, freeParms): - shape is shape of Jacobian [# parms x # free parms] - elts gives the 2D indices of nonzero elements of Jacobian - factors gives the values of the nonzero elements - + """ - + for j, src_j_free in self.jacDoubletPatches: self.jacFactors[j] = freeParms[src_j_free] - + return ((self.nParms, self.nFreeParms), self.jacElts, self.jacFactors) - + @staticmethod @jit(nopython=True, nogil=True) @@ -150,16 +150,16 @@ def _matvec(J_S, v, w): w : :class: `np.ndarray` [# total parameters] (output) """ - + shape, elts, factors = J_S - + for j in range(shape[0]): # total params w[j] = 0. - + for i, (dst, src) in enumerate(elts): w[dst] += factors[i] * v[src] - - + + @staticmethod @jit(nopython=True, nogil=True) def _add_rmatvec(J_S, v, w): @@ -173,11 +173,11 @@ def _add_rmatvec(J_S, v, w): Sparse Jacobian computed by getJacobian() v : :class: `np.ndarray` [# total parameters] w : :class: `np.ndarray` [# free parameters] (output) - + """ _, elts, factors = J_S - + for i, (dst, src) in enumerate(elts): w[src] += factors[i] * v[dst] @@ -213,44 +213,44 @@ def _precomputeMapping(self, isFree, pFree, Source parameters for all doublet ratios. """ - - # by default, assume parameters are fixed + + # by default, assume parameters are fixed sources = np.full(self.nParms, -1, dtype=np.int32) factors = np.zeros(self.nParms, dtype=np.float64) - + # record all free parameters sources[isFree] = pFree[isFree] factors[isFree] = 1. - + for j, src_j in enumerate(tiedSources): if src_j != -1 and isFree[src_j]: # j is tied to src_j with factor tiedFactors[j] sources[j] = pFree[src_j] factors[j] = tiedFactors[j] - + doubletPatches = [] for j, src_j in zip(doubletTargets, doubletSources): if isFree[j] and isFree[src_j]: # j's factor should be v[ p[src_j] ], so that its value # becomes v[ p[j] ] * v[ p[src_j] ]. We will patch it # dynamically at mapping time. - + doubletPatches.append((j, pFree[src_j])) - + self.sources = sources self.factors = factors - + # if there are no patches, we need to create a null array of # the right type and shape to appease Numba if len(doubletPatches) == 0: self.doubletPatches = np.empty((0,2), dtype=np.int32) else: self.doubletPatches = np.array(doubletPatches, dtype=np.int32) - + # record fixed parameter mask self.isFixed = (self.sources == -1) - + def _precomputeJacobian(self): """ Precompute and store as much of the Jacobian of the @@ -262,11 +262,11 @@ def _precomputeJacobian(self): Jacobian. """ - + isLive = np.logical_not(self.isFixed) liveIdx = np.where(isLive)[0] nLive = len(liveIdx) - + # jacElts compresses the source/factor arrays to a sparse array # of coefficients for just the live (non-fixed) # parameters, plus second coeffs for each ratio param of a @@ -276,32 +276,32 @@ def _precomputeJacobian(self): nElts = nLive + nDoublets jacElts = np.empty((nElts,2), dtype=np.int32) jacFactors = np.empty(nElts, dtype=self.factors.dtype) - + jacElts[:nLive,0] = liveIdx jacElts[:nLive,1] = self.sources[liveIdx] jacFactors[:nLive] = self.factors[liveIdx] # makes a copy - + # Create doublet patches for Jacobian w/r to compressed array, # and record new patch locations - + # offsets of orig parms in live array liveOffsets = np.cumsum(isLive) - isLive - + jacDoubletPatches = np.empty((2*len(self.doubletPatches), 2), dtype=np.int32) - + for i, (j, p_src_j) in enumerate(self.doubletPatches): - + p_j = self.sources[j] - + # jacElts already has coefficient (j, p[j]). # its factor should be patched from v[ p[src_j] ] jacDoubletPatches[2*i, :] = (liveOffsets[j], p_src_j) - + # add a second coefficient (j, p[src_j]). # its factor should be patched from v[ p[j] ] jacElts[nLive + i,:] = (j, p_src_j) jacDoubletPatches[2*i+1,:] = (nLive + i, p_j) - + self.jacElts = jacElts self.jacFactors = jacFactors self.jacDoubletPatches = jacDoubletPatches diff --git a/py/fastspecfit/emline_fit/sparse_rep.py b/py/fastspecfit/emline_fit/sparse_rep.py index 55ec92ba..b63a2707 100644 --- a/py/fastspecfit/emline_fit/sparse_rep.py +++ b/py/fastspecfit/emline_fit/sparse_rep.py @@ -19,7 +19,7 @@ class EMLineJacobian(LinearOperator): M is the camera's resolution matrix W is a diagonal matrix of per-observation weights - Note that we precompute W*M*J_I for each camera + Note that we precompute W*M*J_I for each camera with the external mulWMJ() function. This product has one contiguous run of nonzero entries per column, while J_S has either one or two nonzero entries @@ -35,7 +35,7 @@ class EMLineJacobian(LinearOperator): """ def __init__(self, shape, nLineFreeParms, camerapix, jacs, - J_S, J_P=None): + J_S, J_P=None): """ Parameters ---------- @@ -51,21 +51,21 @@ def __init__(self, shape, nLineFreeParms, camerapix, jacs, Parameter expansion Jacobian. J_P : :class:`np.ndarray` [# patch params x # wavelength bins] or :None: Jacobian of patch pedestals, if used. - + """ self.camerapix = camerapix self.jacs = tuple(jacs) self.J_S = J_S self.J_P = J_P self.nLineFreeParms = nLineFreeParms - + if J_P is not None: nPatchParms = J_P[0].shape[0] - + dtype = jacs[0][1].dtype nParms = jacs[0][1].shape[0] # num line params in full set self.vFull = np.empty(nParms, dtype=dtype) - + super().__init__(dtype, shape) @@ -82,18 +82,18 @@ def _matvec(self, v): w : :class:`np.ndarray [# of wavelength bins] """ - + nBins = self.shape[0] w = np.empty(nBins, dtype=v.dtype) vLines = v[:self.nLineFreeParms] - + # return result in self.vFull ParamsMapping._matvec(self.J_S, vLines, self.vFull) - + for campix, jac in zip(self.camerapix, self.jacs): s, e = campix - + # write result to w[s:e] self._matvec_J(jac, self.vFull, w[s:e]) @@ -101,7 +101,7 @@ def _matvec(self, v): if self.J_P is not None: vPatches = v[self.nLineFreeParms:] _matvec_J_add(self.J_P, vPatches, w) - + return w @@ -121,7 +121,7 @@ def _matmat(self, M): w : :class:`np.ndarray [N x # of wavelength bins] """ - + nBins = self.shape[0] nVecs = M.shape[1] R = np.empty((nVecs, nBins), dtype=M.dtype) # transpose of result @@ -130,7 +130,7 @@ def _matmat(self, M): w = R[i,:] v = M[:,i].ravel() - + # return result in self.vFull vLines = v[:self.nLineFreeParms] ParamsMapping._matvec(self.J_S, vLines, self.vFull) @@ -140,15 +140,15 @@ def _matmat(self, M): # write result to w[s:e] self._matvec_J(jac, self.vFull, w[s:e]) - + # add contribution of patch pedestals, if any if self.J_P is not None: vPatches = v[self.nLineFreeParms:] _matvec_J_add(self.J_P, vPatches, w) - + return R.T - - + + def _rmatvec(self, v): """ Compute right matrix-vector product v * J.T, where we are J. @@ -165,23 +165,23 @@ def _rmatvec(self, v): nFreeParms = self.shape[1] w = np.zeros(nFreeParms, dtype=v.dtype) - + wLines = w[:self.nLineFreeParms] for campix, jac in zip(self.camerapix, self.jacs): s, e = campix - + # return result in self.vFull self._rmatvec_J(jac, v[s:e], self.vFull) - + # add result to w ParamsMapping._add_rmatvec(self.J_S, self.vFull, wLines) - + if self.J_P is not None: wPatches = w[self.nLineFreeParms:] self._rmatvec_J(self.J_P, v, wPatches) - + return w - + # # Multiply ideal Jacobian J * v, writing result to w. @@ -198,16 +198,16 @@ def _matvec_J(J, v, w): J : :class:`np.ndarray` [# wavelength bins x # line parameters] v : :class:`np.ndarray` [# of line parameters] w : :class:`np.ndarray` [# wavength bins] (output param) - + """ - + nbins = len(w) for j in range(nbins): w[j] = 0. - + _matvec_J_add(J, v, w) - - + + @staticmethod @jit(nopython=True, nogil=True) def _rmatvec_J(J, v, w): @@ -218,18 +218,18 @@ def _rmatvec_J(J, v, w): Parameters ---------- J : :class:`np.ndarray` [# wavelength bins x # line parameters] - v : :class:`np.ndarray` [# wavength bins] + v : :class:`np.ndarray` [# wavength bins] w : :class:`np.ndarray` [# of line parameters] (output param) """ - + endpts, values = J nvars = endpts.shape[0] - + for i in range(nvars): s, e = endpts[i] vals = values[i] # row i of transpose - + acc = 0. for j in range(e - s): acc += vals[j] * v[j + s] @@ -241,7 +241,7 @@ def _matvec_J_add(J, v, w): """ Multiply partial Jacobian J with v, *adding* result to supplied w. - + Parameters ---------- J : :class:`np.ndarray` [# wavelength bins x # line parameters] @@ -253,15 +253,15 @@ def _matvec_J_add(J, v, w): This function is standalone, rather than a static method of the EMLineJacobian class, only because Numba cannot call a static class method from JIT'd code. - + """ - + endpts, values = J nvars = endpts.shape[0] - + for i in range(nvars): s, e = endpts[i] vals = values[i] # column i - + for j in range(e - s): - w[j + s] += vals[j] * v[i] + w[j + s] += vals[j] * v[i] diff --git a/py/fastspecfit/emline_fit/utils.py b/py/fastspecfit/emline_fit/utils.py index 2813de47..b46d204d 100644 --- a/py/fastspecfit/emline_fit/utils.py +++ b/py/fastspecfit/emline_fit/utils.py @@ -6,7 +6,7 @@ from fastspecfit.util import C_LIGHT -# Do not bother computing normal PDF/CDF if more than this many +# Do not bother computing normal PDF/CDF if more than this many # standard deviations from mean. MAX_SDEV = 5. @@ -16,9 +16,9 @@ def norm_pdf(a): """ PDF of standard normal distribution at a point a """ - + SQRT_2PI = np.sqrt(2 * np.pi) - + return 1/SQRT_2PI * np.exp(-0.5 * a**2) @@ -29,13 +29,13 @@ def norm_cdf(a): """ SQRT1_2 = 1.0 / np.sqrt(2) - + z = np.abs(a) # Optimization (currently disabled because it is not needed): If # |a| > MAX_SDEV, treat the value as extreme and return 0 or 1 as # appropriate. - + #if z > MAX_SDEV: # short-circuit extreme values # if a > 0: # y = 1. @@ -47,7 +47,7 @@ def norm_cdf(a): y = 0.5 * erfc(z * SQRT1_2) if a > 0: y = 1.0 - y - + return y @@ -69,7 +69,7 @@ def max_buffer_width(log_obs_bin_edges, line_sigmas, padding=0): Padding parameter to add to width for future use. """ - + # Find the largest width sigma for any line, and # allocate enough space for twice that much width # in bins, given the smallest observed bin width. From ea8db896c33e822a59b9c1e1bb6acfb90c1faae6 Mon Sep 17 00:00:00 2001 From: jdbuhler Date: Thu, 5 Sep 2024 14:26:23 -0400 Subject: [PATCH 17/20] * remove trailing whitespace from lines in all files * remove shell commands from tops of files that do not contain main or test for __main__ --- py/fastspecfit/continuum.py | 4 +- py/fastspecfit/emlines.py | 12 +++--- py/fastspecfit/fastspecfit.py | 21 ++++++++-- py/fastspecfit/igm.py | 3 +- py/fastspecfit/io.py | 13 +++--- py/fastspecfit/linemasker.py | 5 +-- py/fastspecfit/qa.py | 50 +++++++++++------------ py/fastspecfit/resolution.py | 42 +++++++++---------- py/fastspecfit/test/build_test_spectra.py | 5 +-- py/fastspecfit/util.py | 26 ++++++------ py/fastspecfit/webapp/settings.py | 2 - py/fastspecfit/webapp/urls.py | 2 - 12 files changed, 95 insertions(+), 90 deletions(-) diff --git a/py/fastspecfit/continuum.py b/py/fastspecfit/continuum.py index c3525c41..f42742f5 100644 --- a/py/fastspecfit/continuum.py +++ b/py/fastspecfit/continuum.py @@ -857,13 +857,13 @@ def _stellar_objective(self, params, templateflux, spec_resid = resid[:resid_split] np.subtract(modelflux, specflux, spec_resid) spec_resid *= specistd - + if synthphot: modelflam = self.continuum_to_photometry(fullmodel) phot_resid = resid[resid_split:] np.subtract(modelflam, objflam, phot_resid) phot_resid *= objflamistd - + return resid diff --git a/py/fastspecfit/emlines.py b/py/fastspecfit/emlines.py index 90e7ccfa..e5bafe8e 100644 --- a/py/fastspecfit/emlines.py +++ b/py/fastspecfit/emlines.py @@ -134,7 +134,7 @@ def compute_inrange_lines(self, redshift, wavelims=(3600., 9900.), wavepad=5*0.8 ((zlinewave > (wavelims[0] + wavepad)) & \ (zlinewave < (wavelims[1] - wavepad))) - + def build_linemodels(self, separate_oiii_fit=True): """Build emission line model tables, with and without suppression of broad lines. Establish fixed (i.e., forced to @@ -142,25 +142,25 @@ def build_linemodels(self, separate_oiii_fit=True): relationships for each model, as well as the relationship between lines and their parameters, which we record in the line table. - + Parameter fixing needs to know which lines are within the observed wavelength ranges of the cameras, so we first add this information to the line table. """ - + def create_model(tying_info, forceFixed=[]): """Given the tying info for the model and the list of in-range lines, determine which parameters of the model are fixed and free, and create the model's table. - + We fix all parameters for out-of-range lines to zero, unless the parameter is tied to a parameter for an in-range line. We also allow the caller to specify parameters that should be fixed regardless. """ - + n_params = len(self.param_table) isfixed = np.full(n_params, False, dtype=bool) @@ -355,7 +355,7 @@ def tie_line(tying_info, line_params, source_linename, amp_factor=None): def summarize_linemodel(self, linemodel): """Simple function to summarize an input linemodel.""" - + def _print(line_mask): for line in np.where(line_mask)[0]: line_name = self.line_table['name'][line] diff --git a/py/fastspecfit/fastspecfit.py b/py/fastspecfit/fastspecfit.py index 8509d1b6..1d7bb64b 100644 --- a/py/fastspecfit/fastspecfit.py +++ b/py/fastspecfit/fastspecfit.py @@ -1,4 +1,3 @@ -#!/usr/bin/env python """ fastspecfit.fastspecfit ======================= @@ -15,6 +14,10 @@ from fastspecfit.singlecopy import sc_data from fastspecfit.util import BoxedScalar, MPPool +import cProfile as profile +import pstats +pr = profile.Profile() +shot = 0 def fastspec_one(iobj, data, out_dtype, broadlinefit=True, fastphot=False, constrain_age=False, no_smooth_continuum=False, @@ -41,11 +44,20 @@ def fastspec_one(iobj, data, out_dtype, broadlinefit=True, fastphot=False, # output structure out = BoxedScalar(out_dtype) + global shot + if shot > 0: + pr.enable() + continuummodel, smooth_continuum = continuum_specfit( data, out, templates, igm, phot, constrain_age=constrain_age, no_smooth_continuum=no_smooth_continuum, fastphot=fastphot, debug_plots=debug_plots) + if shot > 0: + pr.disable() + else: + shot = 1 + # Optionally fit the emission-line spectrum. if fastphot: emmodel = None @@ -212,6 +224,10 @@ def fastspec(fastphot=False, stackfit=False, args=None, comm=None, verbose=False log.info(f'Fitting {Spec.ntargets} object(s) took {time.time()-t0:.2f} seconds.') + st = pstats.Stats(pr).strip_dirs().sort_stats("cumulative") + st.print_stats() + st.print_callees() + write_fastspecfit(results, meta, modelspectra=modelspectra, outfile=args.outfile, specprod=Spec.specprod, coadd_type=Spec.coadd_type, fphotofile=sc_data.photometry.fphotofile, @@ -235,7 +251,7 @@ def fastphot(args=None, comm=None): Parameters ---------- args : :class:`argparse.Namespace` or ``None`` - Required and optional arguments parsed via inputs to the command line. + Required and optional arguments parsed via inputs to the command line. comm : :class:`mpi4py.MPI.MPI.COMM_WORLD` or `None` Intracommunicator used with MPI parallelism. @@ -606,4 +622,3 @@ def create_output_table(result_records, meta, units, stackfit=False): results = hstack((results, Table(np.array(result_records), units=units))) return results - diff --git a/py/fastspecfit/igm.py b/py/fastspecfit/igm.py index 3dad887d..60ca7c45 100644 --- a/py/fastspecfit/igm.py +++ b/py/fastspecfit/igm.py @@ -1,4 +1,3 @@ -#!/usr/bin/env python """ fastspecfit.igm =============== @@ -198,7 +197,7 @@ def _tLCDLA(zS, lobs): 0.04696 * z**3 - \ 0.01779 * z**3.3 * (lobs[i]/lamL) ** -0.3 - \ 0.1347 * (lobs[i]/lamL) ** 2 - \ - 0.2905 * (lobs[i]/lamL) ** -0.3 + 0.2905 * (lobs[i]/lamL) ** -0.3 else: r[i] = \ 0.04696 * z**3 - \ diff --git a/py/fastspecfit/io.py b/py/fastspecfit/io.py index 904a284e..ab7b3466 100644 --- a/py/fastspecfit/io.py +++ b/py/fastspecfit/io.py @@ -1,4 +1,3 @@ -#!/usr/bin/env python """ fastspecfit.io ============== @@ -36,8 +35,8 @@ # fibermap and exp_fibermap columns to read FMCOLS = ('TARGETID', 'TARGET_RA', 'TARGET_DEC', 'COADD_FIBERSTATUS', 'OBJTYPE', 'PHOTSYS', 'RELEASE', 'BRICKNAME', 'BRICKID', 'BRICK_OBJID', - #'FIBERFLUX_G', 'FIBERFLUX_R', 'FIBERFLUX_Z', - #'FIBERTOTFLUX_G', 'FIBERTOTFLUX_R', 'FIBERTOTFLUX_Z', + #'FIBERFLUX_G', 'FIBERFLUX_R', 'FIBERFLUX_Z', + #'FIBERTOTFLUX_G', 'FIBERTOTFLUX_R', 'FIBERTOTFLUX_Z', #'FLUX_G', 'FLUX_R', 'FLUX_Z', 'FLUX_W1', 'FLUX_W2', #'FLUX_IVAR_G', 'FLUX_IVAR_R', 'FLUX_IVAR_Z', 'FLUX_IVAR_W1', 'FLUX_IVAR_W2' ) @@ -945,7 +944,7 @@ def one_spectrum(iobj, specdata, meta, ebv, fastphot, specdata['mask'].append(specdata['mask0'][icam]) specdata['res'].append(Resolution(res)) - + if len(cameras) == 0: errmsg = 'No good data, which should never happen.' log.critical(errmsg) @@ -956,7 +955,7 @@ def one_spectrum(iobj, specdata, meta, ebv, fastphot, for key in ('wave', 'flux', 'ivar', 'mask', 'res'): del specdata[key + '0'] specdata[key] = tuple(specdata[key]) - + # Pre-compute some convenience variables for "un-hstacking" # an "hstacked" spectrum. specdata['cameras'] = np.array(cameras) @@ -968,7 +967,7 @@ def one_spectrum(iobj, specdata, meta, ebv, fastphot, specdata['camerapix'] = np.zeros((ncam, 2), np.int16) specdata['camerapix'][:, 0] = c_starts specdata['camerapix'][:, 1] = c_ends - + # use the coadded spectrum to build a robust emission-line mask LM = LineMasker(emline_table) pix = LM.build_linemask( @@ -1351,7 +1350,7 @@ def one_stacked_spectrum(iobj, specdata, synthphot): for key in ('wave', 'flux', 'ivar', 'mask', 'res'): del specdata[key + '0'] specdata[key] = tuple(specdata[key]) - + # Pre-compute some convenience variables for "un-hstacking" # an "hstacked" spectrum. specdata['cameras'] = np.array(cameras) diff --git a/py/fastspecfit/linemasker.py b/py/fastspecfit/linemasker.py index 14bea576..d4ebe5d9 100644 --- a/py/fastspecfit/linemasker.py +++ b/py/fastspecfit/linemasker.py @@ -84,7 +84,7 @@ def _get_contpix(zlinewaves, sigmas, nsigma_factor=2., linemask=None, lya=False) for linename, zlinewave, sigma in zip(linenames, zlinewaves, linesigmas_ang): # skip fixed (e.g., hbeta_broad) lines if sigma <= 0.: - continue + continue I = _get_linepix(zlinewave, sigma) if len(I) > 0: linemask[I] = True @@ -424,7 +424,7 @@ def _fit_patches(continuum_patches, patchMap, linemodel, debug_plots=False, else: isNarrow = EMFit.isNarrow * Ifree if np.any(isNarrow): - linesigma_narrow = np.atleast_1d(linesigmas[isNarrow])[0] + linesigma_narrow = np.atleast_1d(linesigmas[isNarrow])[0] linevshift_narrow = np.atleast_1d(linevshifts[isNarrow])[0] maxsnr_narrow = np.max(linesnrs[isNarrow]) else: @@ -785,4 +785,3 @@ def _niceline(line): } return out - diff --git a/py/fastspecfit/qa.py b/py/fastspecfit/qa.py index afdfaaa1..9fccf0fd 100644 --- a/py/fastspecfit/qa.py +++ b/py/fastspecfit/qa.py @@ -1,4 +1,3 @@ -#!/usr/bin/env python """ fastspecfit.qa ============== @@ -300,7 +299,7 @@ def major_formatter(x, pos): elif fastspec['HGAMMA_BROAD_AMP']*np.sqrt(fastspec['HGAMMA_BROAD_AMP_IVAR']) > snrcut: leg_broad['ewbalmer_broad'] = r'EW(H$\gamma)_{\mathrm{broad}}=$'+r'${:.1f}$'.format(fastspec['HGAMMA_BROAD_EW'])+r' $\AA$' - if (fastspec['OII_3726_AMP']*np.sqrt(fastspec['OII_3726_AMP_IVAR']) > snrcut or + if (fastspec['OII_3726_AMP']*np.sqrt(fastspec['OII_3726_AMP_IVAR']) > snrcut or fastspec['OII_3729_AMP']*np.sqrt(fastspec['OII_3729_AMP_IVAR']) > snrcut): leg_narrow['ewoii'] = r'EW([OII])'+r'$={:.1f}$'.format(fastspec['OII_3726_EW']+fastspec['OII_3729_EW'])+r' $\AA$' @@ -315,24 +314,24 @@ def major_formatter(x, pos): elif fastspec['HGAMMA_AMP']*np.sqrt(fastspec['HGAMMA_AMP_IVAR']) > snrcut: leg_narrow['ewbalmer_narrow'] = r'EW(H$\gamma)=$'+r'${:.1f}$'.format(fastspec['HGAMMA_EW'])+r' $\AA$' - if (fastspec['HALPHA_AMP']*np.sqrt(fastspec['HALPHA_AMP_IVAR']) > snrcut and + if (fastspec['HALPHA_AMP']*np.sqrt(fastspec['HALPHA_AMP_IVAR']) > snrcut and fastspec['HBETA_AMP']*np.sqrt(fastspec['HBETA_AMP_IVAR']) > snrcut): leg_narrow['hahb'] = r'$\mathrm{H}\alpha/\mathrm{H}\beta=$'+r'${:.3f}$'.format(fastspec['HALPHA_FLUX']/fastspec['HBETA_FLUX']) - if 'hahb' not in leg_narrow.keys() and (fastspec['HBETA_AMP']*np.sqrt(fastspec['HBETA_AMP_IVAR']) > snrcut and + if 'hahb' not in leg_narrow.keys() and (fastspec['HBETA_AMP']*np.sqrt(fastspec['HBETA_AMP_IVAR']) > snrcut and fastspec['HGAMMA_AMP']*np.sqrt(fastspec['HGAMMA_AMP_IVAR']) > snrcut): leg_narrow['hbhg'] = r'$\mathrm{H}\beta/\mathrm{H}\gamma=$'+r'${:.3f}$'.format(fastspec['HBETA_FLUX']/fastspec['HGAMMA_FLUX']) - if (fastspec['HBETA_AMP']*np.sqrt(fastspec['HBETA_AMP_IVAR']) > snrcut and - fastspec['OIII_5007_AMP']*np.sqrt(fastspec['OIII_5007_AMP_IVAR']) > snrcut and + if (fastspec['HBETA_AMP']*np.sqrt(fastspec['HBETA_AMP_IVAR']) > snrcut and + fastspec['OIII_5007_AMP']*np.sqrt(fastspec['OIII_5007_AMP_IVAR']) > snrcut and fastspec['HBETA_FLUX'] > 0 and fastspec['OIII_5007_FLUX'] > 0): leg_narrow['oiiihb'] = r'$\log_{10}(\mathrm{[OIII]/H}\beta)=$'+r'${:.3f}$'.format( np.log10(fastspec['OIII_5007_FLUX']/fastspec['HBETA_FLUX'])) - if (fastspec['HALPHA_AMP']*np.sqrt(fastspec['HALPHA_AMP_IVAR']) > snrcut and - fastspec['NII_6584_AMP']*np.sqrt(fastspec['NII_6584_AMP_IVAR']) > snrcut and + if (fastspec['HALPHA_AMP']*np.sqrt(fastspec['HALPHA_AMP_IVAR']) > snrcut and + fastspec['NII_6584_AMP']*np.sqrt(fastspec['NII_6584_AMP_IVAR']) > snrcut and fastspec['HALPHA_FLUX'] > 0 and fastspec['NII_6584_FLUX'] > 0): leg_narrow['niiha'] = r'$\log_{10}(\mathrm{[NII]/H}\alpha)=$'+r'${:.3f}$'.format( np.log10(fastspec['NII_6584_FLUX']/fastspec['HALPHA_FLUX'])) - if (fastspec['OII_3726_AMP']*np.sqrt(fastspec['OII_3726_AMP_IVAR']) > snrcut or + if (fastspec['OII_3726_AMP']*np.sqrt(fastspec['OII_3726_AMP_IVAR']) > snrcut or fastspec['OII_3729_AMP']*np.sqrt(fastspec['OII_3729_AMP_IVAR']) > snrcut): #if fastspec['OII_DOUBLET_RATIO'] != 0: leg_narrow['oii_doublet'] = r'[OII] $\lambda3726/\lambda3729={:.3f}$'.format(fastspec['OII_DOUBLET_RATIO']) @@ -347,7 +346,7 @@ def major_formatter(x, pos): sedmodel, sedphot = CTools.templates2data( templates.flux, templates.wave, redshift=redshift, dluminosity=dlum, photsys=metadata['PHOTSYS'], - synthphot=True, coeff=fastspec['COEFF'] * CTools.massnorm, + synthphot=True, coeff=fastspec['COEFF'] * CTools.massnorm, get_abmag=True) else: sedmodel = CTools.build_stellar_continuum( @@ -398,7 +397,7 @@ def major_formatter(x, pos): else: contmodel = CTools.build_stellar_continuum( templates.flux_nolines, fastspec['COEFF'], - vdisp=fastspec['VDISP'], conv_pre=templates.conv_pre_nolines, + vdisp=fastspec['VDISP'], conv_pre=templates.conv_pre_nolines, ebv=fastspec['AV'] / Templates.klambda(5500.) ) @@ -421,7 +420,7 @@ def major_formatter(x, pos): fullsmoothcontinuum = np.zeros_like(fullwave) else: fullsmoothcontinuum = CTools.smooth_continuum( - fullwave, np.hstack(desiresiduals), np.hstack(data['ivar']), + fullwave, np.hstack(desiresiduals), np.hstack(data['ivar']), np.hstack(data['linemask']), camerapix=data['camerapix']) desismoothcontinuum = [] @@ -460,7 +459,7 @@ def major_formatter(x, pos): if not os.path.isfile(cutoutjpeg): wait = 5 # seconds cmd = 'wget -q -O {outfile} https://www.legacysurvey.org/viewer/jpeg-cutout?ra={ra}&dec={dec}&width={width}&height={height}&layer={layer}' - cmd = cmd.format(outfile=cutoutjpeg, ra=metadata['RA'], dec=metadata['DEC'], + cmd = cmd.format(outfile=cutoutjpeg, ra=metadata['RA'], dec=metadata['DEC'], width=width, height=height, layer=layer) log.info(cmd) try: @@ -646,7 +645,7 @@ def major_formatter(x, pos): sedmodel_abmag = -2.5*np.log10(sedmodel * factor) sedax.plot(sedwave / 1e4, sedmodel_abmag, color='slategrey', alpha=0.9, zorder=1) - sedax.scatter(sedphot['lambda_eff']/1e4, sedphot['abmag'], marker='s', + sedax.scatter(sedphot['lambda_eff']/1e4, sedphot['abmag'], marker='s', s=400, color='k', facecolor='none', linewidth=2, alpha=1.0, zorder=2) if not fastphot: @@ -694,11 +693,11 @@ def major_formatter(x, pos): raise('Problem here!') #print(sed_ymin, sed_ymax) - #sedax.set_xlabel(r'Observed-frame Wavelength ($\mu$m)') + #sedax.set_xlabel(r'Observed-frame Wavelength ($\mu$m)') sedax.set_xlim(phot_wavelims[0], phot_wavelims[1]) sedax.set_xscale('log') - sedax.set_ylabel('AB mag') - #sedax.set_ylabel(r'Apparent Brightness (AB mag)') + sedax.set_ylabel('AB mag') + #sedax.set_ylabel(r'Apparent Brightness (AB mag)') sedax.set_ylim(sed_ymin, sed_ymax) sedax.xaxis.set_major_formatter(major_formatter) @@ -712,7 +711,7 @@ def major_formatter(x, pos): sedax_twin = sedax.twiny() sedax_twin.set_xlim(phot_wavelims[0]/(1+redshift), phot_wavelims[1]/(1+redshift)) sedax_twin.set_xscale('log') - #sedax_twin.set_xlabel(r'Rest-frame Wavelength ($\mu$m)') + #sedax_twin.set_xlabel(r'Rest-frame Wavelength ($\mu$m)') sedax_twin.xaxis.set_major_formatter(major_formatter) restticks = np.array([0.02, 0.05, 0.1, 0.2, 0.5, 1.0, 1.5, 3.0, 5.0, 10.0, 15.0, 20.0]) @@ -865,7 +864,7 @@ def major_formatter(x, pos): thislinename = thislinename.replace('-', ' ') if 'CIII]' in thislinename: thislinename = thislinename+'\n' - specax.text(oneline['restwave']*(1+redshift)/1e4, spec_ymax*0.97, thislinename, + specax.text(oneline['restwave']*(1+redshift)/1e4, spec_ymax*0.97, thislinename, ha='center', va='top', rotation=270, fontsize=12, alpha=0.5) else: if '4363' in linename: @@ -932,7 +931,7 @@ def major_formatter(x, pos): plotline[s:e] = oneline # stupid hack; where cameras overlap (e.g., # Halpha on sv1-bright-22923-39627731570268174), - # the wavelengths are out of order. + # the wavelengths are out of order. srt = np.argsort(fullwave) xx.plot(fullwave[srt] / 1e4, plotline[srt], lw=1, alpha=0.8, color=col2[icam]) @@ -1131,23 +1130,23 @@ def major_formatter(x, pos): _ = leg_broad.pop('dv_broad') ibox = 0 - #fig.text(startpos, toppos, '\n'.join(txt), ha='left', va='top', fontsize=legfntsz, - fig.text(boxpos[ibox], toppos, '\n'.join(txt), ha='left', va='top', fontsize=legfntsz, + #fig.text(startpos, toppos, '\n'.join(txt), ha='left', va='top', fontsize=legfntsz, + fig.text(boxpos[ibox], toppos, '\n'.join(txt), ha='left', va='top', fontsize=legfntsz, bbox=bbox, linespacing=1.4) ibox += 1 txt = [ r'{}'.format(leg['absmag_r']), r'{}'.format(leg['absmag_gr']), - r'{}'.format(leg['absmag_rz']), + r'{}'.format(leg['absmag_rz']), '', r'{}'.format(leg['dn4000_model']) ] if 'dn4000_spec' in legkeys: txt += [r'{}'.format(leg['dn4000_spec'])] - #fig.text(startpos+deltapos, toppos, '\n'.join(txt), ha='left', va='top', - fig.text(boxpos[ibox], toppos, '\n'.join(txt), ha='left', va='top', + #fig.text(startpos+deltapos, toppos, '\n'.join(txt), ha='left', va='top', + fig.text(boxpos[ibox], toppos, '\n'.join(txt), ha='left', va='top', fontsize=legfntsz, bbox=bbox, linespacing=1.4) ibox += 1 @@ -1488,4 +1487,3 @@ def _wrap_qa(redrockfile, indx=None, stackfit=False): # if multiprocessing, clean up workers mp_pool.close() - diff --git a/py/fastspecfit/resolution.py b/py/fastspecfit/resolution.py index a140986d..2fb79624 100644 --- a/py/fastspecfit/resolution.py +++ b/py/fastspecfit/resolution.py @@ -16,7 +16,7 @@ class Resolution(object): fitting code. """ - + def __init__(self, data): """ Parameters @@ -25,17 +25,17 @@ def __init__(self, data): Array of values for diagonals [ndiag//2 .. -ndiag//2]; note that some entries at the ends of rows are ignored. All other diagonals are assumed to be zero. - + """ ndiag, dim = data.shape self.shape = (dim, dim) self.ndiag = ndiag self.data = data - + self.rows = None # compute on first use - + def rowdata(self): """ Compute sparse row matrix from diagonal representation. @@ -44,13 +44,13 @@ def rowdata(self): ------- 2D array of size [dim x ndiag] with nonzero entries for each row of the matrix. - + """ if self.rows is None: self.rows = self._dia_to_rows(self.data) return self.rows - + def dot(self, v, out=None): """ Compute our dot product with vector v, storing @@ -66,50 +66,50 @@ def dot(self, v, out=None): Returns ------- array of length [dim] holding output. - + """ if out is None: out = np.empty(self.shape[0], dtype=v.dtype) return self._matvec(self.data, v, out) - - + + @staticmethod @jit(nopython=True, fastmath=False, nogil=True) def _matvec(D, v, out): """ Compute matrix-vector product of a Resolution matrix A and a vector v. - + A has shape dim x dim, and its nonzero entries are given by a matrix D of size ndiag x dim, providing values on diagonals [ndiag//2 .. -ndiag//2], inclusive. v is an array of length dim. - + Return result in array out of size dim. - + """ out[:] = 0. - + ndiag, dim = D.shape for i in range(ndiag): k = ndiag//2 - i # diagonal offset - + i_start = np.maximum(0, -k) j_start = np.maximum(0, k) j_end = np.minimum(dim + k, dim) - + for n in range(j_end - j_start): out[i_start + n] += D[i, j_start + n] * v[j_start + n] - + return out - - + + @staticmethod @jit(nopython=True, fastmath=False, nogil=True) def _dia_to_rows(D): """ Convert a diagonally sparse matrix M in the form stored by DESI into a sparse row representation. - + Input M is represented as a 2D array D of size ndiag x nrow, whose rows are M's diagonals: M[i,j] = D[ndiag//2 - (j - i), j] @@ -121,7 +121,7 @@ def _dia_to_rows(D): hdiag = ndiag//2 A = np.empty((dim, ndiag), dtype=D.dtype) - + for i in range(dim): # min and max column for row jmin = np.maximum(i - hdiag, 0) @@ -129,5 +129,5 @@ def _dia_to_rows(D): for j in range(jmin, jmax + 1): A[i, j - (i - hdiag)] = D[hdiag + i - j, j] - + return A diff --git a/py/fastspecfit/test/build_test_spectra.py b/py/fastspecfit/test/build_test_spectra.py index d2aed2ad..4ba44655 100644 --- a/py/fastspecfit/test/build_test_spectra.py +++ b/py/fastspecfit/test/build_test_spectra.py @@ -23,7 +23,7 @@ coaddfile = 'coadd-{}-{}-thru{}.fits'.format(petal, tile, night) redrockfile = coaddfile.replace('coadd-', 'redrock-') -#spec, redrock = read_tile_spectra(tile, night, specprod=specprod, targets=targetid, +#spec, redrock = read_tile_spectra(tile, night, specprod=specprod, targets=targetid, # group='cumulative', redrock=True, coadd=True) coaddfile = os.path.join(datadir, coaddfile) @@ -47,7 +47,7 @@ print('Writing {}'.format(os.path.join('data', os.path.basename(redrockfile)))) write_zbest(os.path.join('data', os.path.basename(redrockfile)), - zbest, fibermap, expfibermap, tsnr2, template_version, + zbest, fibermap, expfibermap, tsnr2, template_version, archetype_version, spec_header=spechdr) spec = read_spectra(coaddfile).select(targets=targetid) @@ -75,4 +75,3 @@ #print('Writing {}'.format(outtractorfile)) #os.makedirs(os.path.dirname(outtractorfile), exist_ok=True) #fitsio.write(outtractorfile, tractor.as_array(), header=tractorhdr, clobber=True) - diff --git a/py/fastspecfit/util.py b/py/fastspecfit/util.py index 4aaba632..b64c67c2 100644 --- a/py/fastspecfit/util.py +++ b/py/fastspecfit/util.py @@ -28,7 +28,7 @@ class BoxedScalar(object): field to unbox the scalar. """ - + def __init__(self, dtype): self.value = np.zeros(1, dtype=dtype)[0] @@ -42,7 +42,7 @@ def __setitem__(self, key, v): class MPPool(object): """ A Pool encapsulates paraallel execution with a - multiprocessing.Pool, falling back to sequential + multiprocessing.Pool, falling back to sequential execution in the current process if just one worker is requested. @@ -51,7 +51,7 @@ class MPPool(object): rather than a list of positional arguments. """ - + def __init__(self, nworkers, initializer=None, init_argdict=None): """ create a pool with nworkers workers, using the current @@ -60,7 +60,7 @@ def __init__(self, nworkers, initializer=None, init_argdict=None): init_argdict on startup in each each worker subprocess. """ - + initfunc = None if initializer is None else self.apply_to_dict if nworkers > 1: @@ -82,7 +82,7 @@ def starmap(self, func, argdicts): as a list of keyword argument dictionaries. """ - + # we cannot pickle a local function, so we must pass # both func and the argument dictionary to the subprocess # worker and have it apply one to the other. @@ -100,7 +100,7 @@ def close(self): """ close our multiprocess pool if we created one """ - + if self.pool is not None: self.pool.close() @@ -196,12 +196,12 @@ def mwdust_transmission(ebv, filtername): 'odin-N419': 4.324, 'odin-N501': 3.540, 'odin-N673': 2.438, - 'hsc2017-g': 3.24, - 'hsc2017-r': 2.276, - 'hsc2017-r2': 2.276, - 'hsc2017-i': 1.633, - 'hsc2017-i2': 1.633, - 'hsc2017-z': 1.263, + 'hsc2017-g': 3.24, + 'hsc2017-r': 2.276, + 'hsc2017-r2': 2.276, + 'hsc2017-i': 1.633, + 'hsc2017-i2': 1.633, + 'hsc2017-z': 1.263, 'hsc2017-y': 1.075, 'suprime-IB427': 4.202, 'suprime-IB464': 3.894, @@ -420,7 +420,7 @@ def trapz_rebin(src_x, src_y, bin_centers, out=None, pre=None): out: :class:`numpy.ndarray` or None if not None, an array of correct size and type that will receive the result of the computation. - If None, a new array will be allocated. + If None, a new array will be allocated. pre : preprocessing data computed by trapz_rebin_pre(), if available. If not None, it must correspond diff --git a/py/fastspecfit/webapp/settings.py b/py/fastspecfit/webapp/settings.py index 3e2af1fa..7fbb6373 100644 --- a/py/fastspecfit/webapp/settings.py +++ b/py/fastspecfit/webapp/settings.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - """Django settings for FastSpecFit project. Generated by 'django-admin startproject' using Django 2.0.6. diff --git a/py/fastspecfit/webapp/urls.py b/py/fastspecfit/webapp/urls.py index 6f56f012..9801e5b1 100644 --- a/py/fastspecfit/webapp/urls.py +++ b/py/fastspecfit/webapp/urls.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - """URL Configuration The `urlpatterns` list routes URLs to views. For more information please see: From 86414a857e20f885c1dfdead4490960927b5c5f7 Mon Sep 17 00:00:00 2001 From: jdbuhler Date: Thu, 5 Sep 2024 14:33:09 -0400 Subject: [PATCH 18/20] remove header comments with license and charset info --- py/fastspecfit/__init__.py | 2 -- py/fastspecfit/test/test_top_level.py | 2 -- 2 files changed, 4 deletions(-) diff --git a/py/fastspecfit/__init__.py b/py/fastspecfit/__init__.py index f11b8e1b..6d95b0b6 100644 --- a/py/fastspecfit/__init__.py +++ b/py/fastspecfit/__init__.py @@ -1,5 +1,3 @@ -# Licensed under a 3-clause BSD style license - see LICENSE.rst -# -*- coding: utf-8 -*- """ fastspecfit =========== diff --git a/py/fastspecfit/test/test_top_level.py b/py/fastspecfit/test/test_top_level.py index f7a79905..2dcd1e31 100644 --- a/py/fastspecfit/test/test_top_level.py +++ b/py/fastspecfit/test/test_top_level.py @@ -1,5 +1,3 @@ -# Licensed under a 3-clause BSD style license - see LICENSE.rst -# -*- coding: utf-8 -*- """ test top-level fastspecfit functions From 0332207bf60b952f6fb3d66c575d29a4dfa333a0 Mon Sep 17 00:00:00 2001 From: jdbuhler Date: Thu, 5 Sep 2024 14:38:22 -0400 Subject: [PATCH 19/20] revert inclusion of profiling code in fastspecfit.py --- py/fastspecfit/fastspecfit.py | 17 ----------------- 1 file changed, 17 deletions(-) diff --git a/py/fastspecfit/fastspecfit.py b/py/fastspecfit/fastspecfit.py index 1d7bb64b..ee545d7c 100644 --- a/py/fastspecfit/fastspecfit.py +++ b/py/fastspecfit/fastspecfit.py @@ -14,10 +14,6 @@ from fastspecfit.singlecopy import sc_data from fastspecfit.util import BoxedScalar, MPPool -import cProfile as profile -import pstats -pr = profile.Profile() -shot = 0 def fastspec_one(iobj, data, out_dtype, broadlinefit=True, fastphot=False, constrain_age=False, no_smooth_continuum=False, @@ -44,20 +40,11 @@ def fastspec_one(iobj, data, out_dtype, broadlinefit=True, fastphot=False, # output structure out = BoxedScalar(out_dtype) - global shot - if shot > 0: - pr.enable() - continuummodel, smooth_continuum = continuum_specfit( data, out, templates, igm, phot, constrain_age=constrain_age, no_smooth_continuum=no_smooth_continuum, fastphot=fastphot, debug_plots=debug_plots) - if shot > 0: - pr.disable() - else: - shot = 1 - # Optionally fit the emission-line spectrum. if fastphot: emmodel = None @@ -224,10 +211,6 @@ def fastspec(fastphot=False, stackfit=False, args=None, comm=None, verbose=False log.info(f'Fitting {Spec.ntargets} object(s) took {time.time()-t0:.2f} seconds.') - st = pstats.Stats(pr).strip_dirs().sort_stats("cumulative") - st.print_stats() - st.print_callees() - write_fastspecfit(results, meta, modelspectra=modelspectra, outfile=args.outfile, specprod=Spec.specprod, coadd_type=Spec.coadd_type, fphotofile=sc_data.photometry.fphotofile, From 78916bc74db79ea38ec4fbaaee50a51658fe2ab9 Mon Sep 17 00:00:00 2001 From: jdbuhler Date: Thu, 5 Sep 2024 14:45:45 -0400 Subject: [PATCH 20/20] extend trailing whitespace and gratuitous header removal to webapp/fastmodel --- py/fastspecfit/webapp/fastmodel/filters.py | 4 +-- py/fastspecfit/webapp/fastmodel/models.py | 24 +++++++------- .../fastmodel/templatetags/my_templatetag.py | 22 ++++++------- py/fastspecfit/webapp/fastmodel/views.py | 33 +++++++++---------- 4 files changed, 39 insertions(+), 44 deletions(-) diff --git a/py/fastspecfit/webapp/fastmodel/filters.py b/py/fastspecfit/webapp/fastmodel/filters.py index 837f433e..6f200850 100644 --- a/py/fastspecfit/webapp/fastmodel/filters.py +++ b/py/fastspecfit/webapp/fastmodel/filters.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - """Custom filters for the FastModel model, which work by selecting FastModel objects in the database based on meeting the desired criteria. @@ -16,7 +14,7 @@ class FastModelFilter(django_filters.FilterSet): """ #field_name is the FastModel object variable - #lookup_expr is used to get ranges (currently using greater/less than or equal to + #lookup_expr is used to get ranges (currently using greater/less than or equal to survey__match = django_filters.CharFilter(field_name='survey', lookup_expr='icontains') program__match = django_filters.CharFilter(field_name='program', lookup_expr='icontains') diff --git a/py/fastspecfit/webapp/fastmodel/models.py b/py/fastspecfit/webapp/fastmodel/models.py index 380a6cd9..f9483b39 100644 --- a/py/fastspecfit/webapp/fastmodel/models.py +++ b/py/fastspecfit/webapp/fastmodel/models.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python - """Each model will be written as a class here, instantiated and populated by load.py, with each model stored as a table in the database and the fields stored as columns. @@ -18,7 +16,7 @@ class FastModel(Model): """ #def __str__(self): # return 'Sample '+self.target_name - + # in FITS table row_index = BigIntegerField(default=0, db_index=True) @@ -124,11 +122,11 @@ class FastModel(Model): smoothcorr_b = FloatField(null=True) smoothcorr_r = FloatField(null=True) smoothcorr_z = FloatField(null=True) - + vdisp = FloatField(null=True) vdisp_ivar = FloatField(null=True) vdisp_err = CharField(max_length=50, default='') - + age = FloatField(null=True) zzsun = FloatField(null=True) logmstar = FloatField(null=True) @@ -160,7 +158,7 @@ class FastModel(Model): absmag10_decam_g = FloatField(null=True) absmag10_decam_r = FloatField(null=True) absmag10_decam_z = FloatField(null=True) - + kcorr00_u = FloatField(null=True) kcorr00_b = FloatField(null=True) kcorr00_v = FloatField(null=True) @@ -246,19 +244,19 @@ class FastModel(Model): narrow_dv = FloatField(null=True) broad_dv = FloatField(null=True) uv_dv = FloatField(null=True) - narrow_dv_err = CharField(max_length=50, default='') - broad_dv_err = CharField(max_length=50, default='') - uv_dv_err = CharField(max_length=50, default='') - + narrow_dv_err = CharField(max_length=50, default='') + broad_dv_err = CharField(max_length=50, default='') + uv_dv_err = CharField(max_length=50, default='') + narrow_sigma = FloatField(null=True) broad_sigma = FloatField(null=True) uv_sigma = FloatField(null=True) narrow_sigmarms = FloatField(null=True) broad_sigmarms = FloatField(null=True) uv_sigmarms = FloatField(null=True) - narrow_sigma_err = CharField(max_length=50, default='') - broad_sigma_err = CharField(max_length=50, default='') - uv_sigma_err = CharField(max_length=50, default='') + narrow_sigma_err = CharField(max_length=50, default='') + broad_sigma_err = CharField(max_length=50, default='') + uv_sigma_err = CharField(max_length=50, default='') mgii_doublet_ratio = FloatField(null=True) oii_doublet_ratio = FloatField(null=True) diff --git a/py/fastspecfit/webapp/fastmodel/templatetags/my_templatetag.py b/py/fastspecfit/webapp/fastmodel/templatetags/my_templatetag.py index 83c91cc3..fefc0a50 100644 --- a/py/fastspecfit/webapp/fastmodel/templatetags/my_templatetag.py +++ b/py/fastspecfit/webapp/fastmodel/templatetags/my_templatetag.py @@ -18,12 +18,12 @@ @register.simple_tag def url_replace(req, field, value): """Replace the old GET value of desired field with a new value. - + Args: req: the http request field: the field to replace value: the new value - + Returns: The updated url with the new value @@ -36,11 +36,11 @@ def url_replace(req, field, value): def url_replace_sort(req, new_sort): """Replace the old GET value of sort with a new value, or negate it if they are equal to sort the opposite way. - + Args: req: the http request new_sort: the sort value a user clicked on - + Returns: The updated url with the new sort value @@ -55,14 +55,14 @@ def url_replace_sort(req, new_sort): else: dict_['sort'] = new_sort return dict_.urlencode() - + @register.simple_tag def url_pull(req): """Return a string describing the search criteria used. - + Args: req: the http request - + Returns: Description of search criteria @@ -80,7 +80,7 @@ def url_pull(req): if dict_["ra__lte"] == "": search += " RA high: 360 |" else: - search += " RA high: " + dict_["ra__lte"] + " |" + search += " RA high: " + dict_["ra__lte"] + " |" entry = True if "dec__gte" in dict_: if dict_["dec__gte"] == "": @@ -116,12 +116,12 @@ def url_pull(req): def viewer_link(ra, dec): """ Creates a string with the viewer link for desired galaxy. - + Args: ra: the ra value to use in link dec: the dec value to use in link - - Returns: + + Returns: Viewer link based on ra and dec values """ baseurl = 'http://legacysurvey.org/viewer/' diff --git a/py/fastspecfit/webapp/fastmodel/views.py b/py/fastspecfit/webapp/fastmodel/views.py index 33ddaa4c..6460c15a 100644 --- a/py/fastspecfit/webapp/fastmodel/views.py +++ b/py/fastspecfit/webapp/fastmodel/views.py @@ -32,11 +32,11 @@ def explore(req): """Returns the explore.html file, or renders the explore.html page after it applies the filter, stores result to session, and sets up pagination. - + Args: req: the http request - - Returns: + + Returns: File stream if user clicked download, otherwise render for explore.html """ @@ -46,7 +46,7 @@ def explore(req): # If download button was pressed return the selected subset of the FITS table. if req.method == 'POST': from fastspecfit.webapp.load import DATADIR, fastspecfile - + print('download: req.GET:', req.GET) query = pickle.loads(req.session['fastmodel_query']) print('Query:', query) @@ -110,7 +110,7 @@ def explore(req): cone_ra = req.GET.get('conera','') cone_dec = req.GET.get('conedec','') cone_rad = req.GET.get('coneradius','') - + # save for form default cone_rad_arcmin = cone_rad if len(cone_ra) and len(cone_dec) and len(cone_rad): @@ -149,12 +149,12 @@ def explore(req): #use pickle to serialize queryset (just the SQL query), and store in session req.session['fastmodel_query'] = pickle.dumps(fastmodel_filtered.query) - + #use django pagination functionality paginator = Paginator(fastmodel_filtered, 50) page_num = req.GET.get('page') page = paginator.get_page(page_num) - + # Include pagination values we will use in html page in the return # statement. for sam in page: @@ -176,7 +176,7 @@ def target(req, target_name): # grab this one (unique) target target = FastModel.objects.all().filter(target_name=target_name) target.order_by('targetid') - + result_index = req.GET.get('index', '-1') try: result_index = int(result_index, 10) @@ -189,7 +189,7 @@ def target(req, target_name): i_prev,_ = get_next_target(req, result_index, direction=-1) has_next = i_next is not None has_prev = i_prev is not None - + return render(req, 'target.html', {'target': target[0], 'target_name': target_name, 'result_index': result_index, @@ -236,13 +236,13 @@ def target_next(req, index): def index(req): """ Renders the homepage from index.html - + Args: req: the http request - - Returns: + + Returns: Render for index.html - + """ return render(req, 'index.html') @@ -303,7 +303,7 @@ def fastmodel_near_radec(ra, dec, rad, tablename='fastmodel', ('SELECT *, ((ux-(%g))*(ux-(%g))+(uy-(%g))*(uy-(%g))+(uz-(%g))*(uz-(%g))) as r2' + ' FROM %s where r2 <= %g %s ORDER BY r2') % (x,x,y,y,z,z, tablename, r2, extra_where)) - + return sample def index(req, **kwargs): @@ -322,7 +322,7 @@ def upload_catalog(req): from django.urls import reverse #from map.views import index from fastspecfit.webapp import settings - + if req.method != 'POST': return HttpResponse('POST only') print('Files:', req.FILES) @@ -344,7 +344,7 @@ def upload_catalog(req): print('Wrote', tmpfn) errtxt = ('%s

Custom catalogs must be a binary FITS binary table with mandatory columns "SURVEY", "PROGRAM", "HEALPIX" and "TARGETID".

') - + #+'See Tips & Tricks for some hints on how to produce such a catalog.

') T = None @@ -426,4 +426,3 @@ def main(): print('Time:', times[i], 'SQL:', queries[i]) #main() -