Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

New resolution matrix implementation #181

Merged
merged 20 commits into from
Sep 16, 2024
Merged

New resolution matrix implementation #181

merged 20 commits into from
Sep 16, 2024

Conversation

jdbuhler
Copy link
Collaborator

@jdbuhler jdbuhler commented Sep 2, 2024

Reimplement enough of desispec.resolution.Resolution to perform matrix-vector operations, which is all we do with it, with a significant speedup via Numba compared to the scipy.sparse implementation.

Retire the old res_emline version of the matrix, which is slower for dot product than either DESI or our new implementation. Allow conversion on the fly from our new impl to the internal rep used by res_emline, which is used way down in the guts of emline fitting but not otherwise. Switching away from res_emline for dot products does result in a small change in the output of emline fitting, which propagates to masking.

In addition, make two small but meaningful changes to least_squares call in continuum fitting: switch from lsmr to exact method for TRF algorithm, which reduces least_squares overhead, and use a looser ftol of 1e-3 instead of 1e-5, which greatly reduces the number of lsq iterations. (These should only be accepted if the output still looks OK!)

In passing, implement several small code cleanups that do not impact the results at all, including more efficient storage of some fields in the data dictionary and avoiding a copy in the _stellar_objective().

jdbuhler and others added 5 commits August 30, 2024 23:19
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.
…ices

* 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
@moustakas
Copy link
Member

Thanks @jdbuhler. The fitting code looks good but running fastqa crashes with

[snip]
Traceback (most recent call last):
  File "/global/homes/i/ioannis/code/desihub/fastspecfit/bin/fastqa", line 9, in <module>
    fastqa(comm=None)
  File "/global/homes/i/ioannis/code/desihub/fastspecfit/py/fastspecfit/qa.py", line 1431, in fastqa
    _wrap_qa(redrockfile, indx)
  File "/global/homes/i/ioannis/code/desihub/fastspecfit/py/fastspecfit/qa.py", line 1407, in _wrap_qa
    for _ in mp_pool.starmap(desiqa_one, qaargs): pass
  File "/global/homes/i/ioannis/code/desihub/fastspecfit/py/fastspecfit/util.py", line 97, in apply_to_dict
    return f(**argdict)
  File "/global/homes/i/ioannis/code/desihub/fastspecfit/py/fastspecfit/qa.py", line 25, in desiqa_one
    qa_fastspec(data, sc_data.templates, fastfit, metadata,
  File "/global/homes/i/ioannis/code/desihub/fastspecfit/py/fastspecfit/qa.py", line 432, in qa_fastspec
    _desiemlines = EMFit.emlinemodel_bestfit(fastspec, fastspec['Z'], np.hstack(data['wave']), data['res_emline'],
KeyError: 'res_emline'
[snip]

Can you take a look and push a fix?

Also, can you please update the doc/changes.rst file with info regarding this PR?

jdbuhler added 7 commits September 3, 2024 10:36
* update documentation of sparse emline Jacobian to std
further limit ftol to 1e-3 for greater speed.
* attempt to clean up some formatting warnings
* add new-style docstrings to emline_fit/utils.py
of functions inner where they have only one caller
* consistently import jit from numba instead of saying '@numba.jit'
@@ -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):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Didn't we want to change fastmath=False here and elsewhere?

@jdbuhler
Copy link
Collaborator Author

jdbuhler commented Sep 5, 2024 via email

Copy link
Member

@moustakas moustakas left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The code changes look fine, but either this PR or #180 has led to non-negligible changes in the continuum-fitting results (as you warned me, off-list).

I added just a couple comments. After you make those changes, let's go ahead and merge this PR and then I'm going to test the fitting results on a larger sample.

@@ -11,7 +11,7 @@


class Inoue14(object):
r"""
"""
IGM absorption from Inoue et al. (2014)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Without the "raw string" formatting, integration tests give a DeprecationWarning. Can we please revert this change?

/global/u2/i/ioannis/code/desihub/fastspecfit/py/fastspecfit/igm.py:14: DeprecationWarning: invalid escape sequence '\m'

modelflux = np.empty(self.wavelen)

for icam, (s, e) in enumerate(camerapix):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's a little (a lot?) OCD, but can we remove the added blank 8 spaces on lines 754 and 756? (My configuration of spacemacs shows these blank characters, so I've been cleaning them up.)

@@ -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
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove blank spaces.

self.convolve = sc_sig.convolve
log.info('Using mkl_fft library for FFTs')
else:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe this should be a debugging comment?

@jdbuhler
Copy link
Collaborator Author

jdbuhler commented Sep 5, 2024 via email

@moustakas
Copy link
Member

Thanks, @jdbuhler, I'm going to go ahead and merge this and then run on a few 10^4 spectra.

@moustakas moustakas merged commit b71cec6 into main Sep 16, 2024
12 checks passed
@moustakas moustakas deleted the new-resmul branch September 16, 2024 17:26
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants