Merge pull request #181 from desihub/new-resmul

New resolution matrix implementation

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

py/fastspecfit/__init__.py

@@ -1,5 +1,3 @@
-# Licensed under a 3-clause BSD style license - see LICENSE.rst
-# -*- coding: utf-8 -*-
 """
 fastspecfit
 ===========

py/fastspecfit/continuum.py

@@ -7,14 +7,14 @@
 """
 import time
 import numpy as np
-import numba
+from numba import jit
 
 from fastspecfit.logger import log
 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):
@@ -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),
@@ -71,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)
 
 
@@ -341,7 +338,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
@@ -592,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,
@@ -631,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,
@@ -748,7 +745,7 @@ 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']
@@ -762,7 +759,7 @@ def continuum_to_spectroscopy(self, contmodel):
  contmodel,
  specwave[icam],
  pre=self.spec_pre[icam])
- modelflux[s:e] = specres[icam] @ resampflux
+ specres[icam].dot(resampflux, out=modelflux[s:e])
 
  return modelflux
 
@@ -834,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,
@@ -860,16 +855,14 @@ 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
 
@@ -990,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
 

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