Merge pull request #47 from lsst/tickets/DM-40988

DM-40988: Resync with upstream

.github/workflows/build.yaml

@@ -28,6 +28,7 @@ jobs:
         run: |
           mamba install -y  "numpy>1.15" scipy "matplotlib>3.1" pandas llvmlite numba "astropy>=3.2" "photutils>=1.7" astroquery coloredlogs scikit-image>=0.20 h5py emcee tqdm mpi4py schwimmbad "iminuit>=2" "coverage>=3.6" configparser coveralls deprecated pyyaml nose rubin-libradtran "getCalspec>=2.0.0"
           # python -c "from getCalspec.rebuild import rebuild_tables; rebuild_tables()"
+          pip install lsst.utils
           mamba install astrometry
 
       - name: Download test data

requirements.txt

@@ -13,11 +13,11 @@ h5py #==2.8.0
 emcee #==3.0rc2
 tqdm
 schwimmbad
-# iminuit>=2
 coverage>=3.6 # <5
 configparser
 coveralls
 deprecated
 pyyaml
 nose
 getCalspec>=2.0.0
+lsst_utils

spectractor/_version.py

@@ -28,5 +28,5 @@
 # OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 # OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
-__version__ = '3.0'
+__version__ = '3.0.1'
 __version_info__ = tuple(map(int, __version__.split('.')))

spectractor/extractor/extractor.py

@@ -100,9 +100,10 @@ def __init__(self, spectrum, amplitude_priors_method="noprior", verbose=False, p
         # This set of fixed parameters was determined so that the reconstructed spectrum has a ZERO bias
         # with respect to the true spectrum injected in the simulation
         # A2 is free only if spectrogram is a simulation or if the order 2/1 ratio is not known and flat
-        params.fixed[params.get_index("A1")] = True  # A1
-        params.fixed[params.get_index("A2")] = (not spectrum.disperser.flat_ratio_order_2over1) and (not ("A2_T" in spectrum.header))
-        params.fixed[params.get_index("A3")] = True  # A3
+        for order in self.diffraction_orders:
+            params.fixed[params.get_index(f"A{order}")] = True
+        if "A2" in params.labels:
+            params.fixed[params.get_index("A2")] = (not spectrum.disperser.flat_ratio_order_2over1) and (not ("A2_T" in spectrum.header))
         params.fixed[params.get_index("D_CCD [mm]")] = True  # D2CCD: spectrogram can not tell something on this parameter: rely on calibrate_spectrum
         params.fixed[params.get_index("shift_x [pix]")] = True  # delta x: if False, extracted spectrum is biased compared with truth
         params.fixed[params.get_index("shift_y [pix]")] = True  # delta y
@@ -199,6 +200,7 @@ def __init__(self, spectrum, amplitude_priors_method="noprior", verbose=False, p
             self.reg = float(spectrum.header['PSF_REG'])
         if self.reg < 0:
             self.reg = parameters.PSF_FIT_REG_PARAM
+        self.trace_r = self.Nx / np.min(self.fwhm_priors)  # spectrophotometric uncertainty principle
         self.my_logger.info(f"\n\tFull forward model fitting with regularisation parameter r={self.reg}.")
         self.Q = np.zeros((self.Nx, self.Nx), dtype="float32")
         self.Q_dot_A0 = np.zeros(self.Nx, dtype="float32")
@@ -278,15 +280,15 @@ def set_mask(self, params=None, fwhmx_clip=3*parameters.PSF_FWHM_CLIP, fwhmy_cli
             profile_params[:, 0] = 1
             profile_params[:, 1] = dispersion_law.real + self.spectrum.spectrogram_x0
             profile_params[:, 2] += dispersion_law.imag - self.bgd_width
-            psf_cube = self.spectrum.chromatic_psf.build_psf_cube(self.pixels, profile_params,
-                                                                  fwhmx_clip=fwhmx_clip,
-                                                                  fwhmy_clip=fwhmy_clip, dtype="float32")
+            psf_cube_masked = self.spectrum.chromatic_psf.build_psf_cube_masked(self.pixels, profile_params,
+                                                                                fwhmx_clip=fwhmx_clip,
+                                                                                fwhmy_clip=fwhmy_clip)
 
-            self.psf_cubes_masked[order] = self.spectrum.chromatic_psf.get_psf_cube_masked(psf_cube, convolve=True)
+            self.psf_cubes_masked[order] = self.spectrum.chromatic_psf.convolve_psf_cube_masked(psf_cube_masked)
             # make rectangular mask per wavelength
             self.boundaries[order], self.psf_cubes_masked[order] = self.spectrum.chromatic_psf.get_boundaries(self.psf_cubes_masked[order])
             self.psf_cube_sparse_indices[order], self.M_sparse_indices[order] = self.spectrum.chromatic_psf.get_sparse_indices(self.psf_cubes_masked[order])
-        mask = np.sum(self.psf_cubes_masked[self.diffraction_orders[0]].reshape(psf_cube.shape[0], psf_cube[0].size), axis=0) == 0
+        mask = np.sum(self.psf_cubes_masked[self.diffraction_orders[0]].reshape(psf_cube_masked.shape[0], psf_cube_masked[0].size), axis=0) == 0
         self.W = np.copy(self.W_before_mask)
         self.W[mask] = 0
         self.sqrtW = sparse.diags(np.sqrt(self.W), format="dia", dtype="float32")
@@ -846,9 +848,10 @@ def run_ffm_minimisation(w, method="newton", niter=2):
         if w.amplitude_priors_method == "spectrum" and w.reg == parameters.PSF_FIT_REG_PARAM:  # pragma: no cover
             my_logger.info("\n\tStart regularization parameter estimation...")
             w_reg = RegFitWorkspace(w, opt_reg=parameters.PSF_FIT_REG_PARAM, verbose=True)
-            w_reg.run_regularisation()
+            w_reg.run_regularisation(Ndof=w.trace_r)
             w.opt_reg = w_reg.opt_reg
             w.reg = np.copy(w_reg.opt_reg)
+            w.trace_r = np.trace(w_reg.resolution)
             w.simulate(*w.params.values)
             if np.trace(w.amplitude_cov_matrix) < np.trace(w.amplitude_priors_cov_matrix):
                 w.my_logger.warning(
@@ -863,7 +866,7 @@ def run_ffm_minimisation(w, method="newton", niter=2):
                 w.amplitude_params_err = np.array(
                     [np.sqrt(w.amplitude_cov_matrix[x, x]) for x in range(w.Nx)])
             w.spectrum.header['PSF_REG'] = w.opt_reg
-            w.spectrum.header['TRACE_R'] = np.trace(w_reg.resolution)
+            w.spectrum.header['TRACE_R'] = w.trace_r
 
         if parameters.DEBUG and parameters.DISPLAY:
             w.plot_fit()
@@ -1322,8 +1325,7 @@ def bgd_model_func(x, y):
     method = "noprior"
     mode = "1D"
 
-    my_logger.info('\n\t  ======================= ChromaticPSF polynomial fit  =============================')
-
+    my_logger.info('\n\t  ======================= ChromaticPSF1D polynomial fit  =============================')
     my_logger.info(f'\n\tStart ChromaticPSF polynomial fit with '
                    f'mode={mode} and amplitude_priors_method={method}...')
     w = s.fit_chromatic_psf(data, bgd_model_func=bgd_model_func, data_errors=err,
@@ -1525,13 +1527,14 @@ def run_spectrogram_deconvolution_psf2d(spectrum, bgd_model_func):
                    f'mode={mode} and amplitude_priors_method={method}...')
     data = spectrum.spectrogram
     err = spectrum.spectrogram_err
+
+    my_logger.info('\n\t  ======================= ChromaticPSF2D polynomial fit  =============================')
     w = s.fit_chromatic_psf(data, bgd_model_func=bgd_model_func, data_errors=err, live_fit=False,
                             amplitude_priors_method=method, mode=mode, verbose=parameters.VERBOSE, analytical=True)
 
     # save results
-    spectrum.spectrogram_fit = s.evaluate(s.set_pixels(mode=mode), poly_params=s.params.values)
-    spectrum.spectrogram_residuals = (data - spectrum.spectrogram_fit - bgd_model_func(np.arange(Nx),
-                                                                                       np.arange(Ny))) / err
+    spectrum.spectrogram_fit = w.model.reshape((w.Ny, w.Nx))  #s.evaluate(s.set_pixels(mode=mode), poly_params=s.params.values)
+    spectrum.spectrogram_residuals = (w.data.reshape((w.Ny, w.Nx)) - spectrum.spectrogram_fit) / w.err.reshape((w.Ny, w.Nx))
     lambdas = spectrum.disperser.grating_pixel_to_lambda(s.get_algebraic_distance_along_dispersion_axis(),
                                                          x0=spectrum.x0, order=spectrum.order)
     s.table['lambdas'] = lambdas

spectractor/extractor/images.py

@@ -1127,11 +1127,13 @@ def find_target_Moffat2D(image, sub_image_subtracted, sub_errors=None):
     XX = np.arange(NX)
     YY = np.arange(NY)
     # find a first guess of the target position
-    avX, sigX = weighted_avg_and_std(XX, np.sum(sub_image_subtracted, axis=0) ** 4)
-    avY, sigY = weighted_avg_and_std(YY, np.sum(sub_image_subtracted, axis=1) ** 4)
+    xprofile = np.sum(sub_image_subtracted, axis=0)
+    yprofile = np.sum(sub_image_subtracted, axis=1)
+    _, sigX = weighted_avg_and_std(XX, xprofile)
+    _, sigY = weighted_avg_and_std(YY, yprofile)
+    avX = np.average(XX, weights=xprofile ** 4)
+    avY = np.average(YY, weights=yprofile ** 4)
     # fit a 2D star profile close to this position
-    # guess = [np.max(sub_image_subtracted),avX,avY,1,1] #for Moffat2Ds
-    # guess = [np.max(sub_image_subtracted),avX-2,avY-2,2,2,0] #for Gauss2D
     psf = Moffat(clip=True)
     total_flux = np.sum(sub_image_subtracted)
     psf.params.values[:3] = [total_flux, avX, avY]
@@ -1141,18 +1143,11 @@ def find_target_Moffat2D(image, sub_image_subtracted, sub_errors=None):
         psf.params.values[1] = avX
         psf.params.values[2] = avY
     mean_prior = 10  # in pixels
-    # bounds = [ [0.5*np.max(sub_image_subtracted),avX-mean_prior,avY-mean_prior,0,-np.inf],
-    # [2*np.max(sub_image_subtracted),avX+mean_prior,avY+mean_prior,np.inf,np.inf] ] #for Moffat2D
-    # bounds = [ [0.5*np.max(sub_image_subtracted),avX-mean_prior,avY-mean_prior,0,0,0],
-    # [100*image.saturation,avX+mean_prior,avY+mean_prior,10,10,np.pi] ] #for Gauss2D
-    # bounds = [[0.5 * np.max(sub_image_subtracted), avX - mean_prior, avY - mean_prior, 2, 0.9 * image.saturation],
-    # [10 * image.saturation, avX + mean_prior, avY + mean_prior, 15, 1.1 * image.saturation]]
-    psf.params.bounds[:3] = [[0.1 * total_flux, 4 * total_flux],
+    psf.params.bounds[:4] = [[0.1 * total_flux, 10 * total_flux],
                              [avX - mean_prior, avX + mean_prior],
-                             [avY - mean_prior, avY + mean_prior]]
+                             [avY - mean_prior, avY + mean_prior],
+                             [0.5*min(sigX, sigY), min(NX, NY)]]
     # fit
-    # star2D = fit_PSF2D(X, Y, sub_image_subtracted, guess=guess, bounds=bounds)
-    # star2D = fit_PSF2D_minuit(X, Y, sub_image_subtracted, guess=guess, bounds=bounds)
     psf.fit_psf(sub_image_subtracted, data_errors=sub_errors, bgd_model_func=image.target_bkgd2D)
     new_avX = psf.params.values[1]
     new_avY = psf.params.values[2]

spectractor/extractor/psf.py

@@ -418,7 +418,8 @@ def evaluate_moffatgauss1d_jacobian(y, amplitude, y_c, gamma, alpha, eta_gauss,
     return J
 
 
-@njit(["float32[:,:](int64[:,:], int64[:,:], float32, float32, float32, float32, float32)"], fastmath=True, cache=True)
+@njit(["float32[:,:](int64[:,:], int64[:,:], float32, float32, float32, float32, float32)",
+       "float64[:,:](float64[:,:], float64[:,:], float32, float32, float32, float32, float32)"], fastmath=True, cache=True)
 def evaluate_moffat2d(x, y, amplitude, x_c, y_c, gamma, alpha):  # pragma: no cover
     r"""Compute a 2D Moffat function, whose integral is normalised to unity.
 

spectractor/extractor/spectrum.py

@@ -28,6 +28,9 @@
 from spectractor.simulation.throughput import TelescopeTransmission
 from spectractor.fit.fitter import FitWorkspace, FitParameters, run_minimisation
 
+from lsst.utils.threads import disable_implicit_threading
+disable_implicit_threading()
+
 
 fits_mappings = {'config': 'CONFIG',
                  'date_obs': 'DATE-OBS',

spectractor/extractor/targets.py

@@ -376,8 +376,9 @@ def load_spectra(self):
                              f"\n\tCalspec? {is_calspec}"
                              f"\n\tNumber of spectra: {len(self.spectra)}"
                              f"\n\tRedshift: {self.redshift}"
-                             f"\n\tEmission spectrum ? {self.emission_spectrum}"
-                             f"\n\tLines: {[l.label for l in self.lines.lines]}")
+                             f"\n\tEmission spectrum ? {self.emission_spectrum}")
+        if self.lines is not None and len(self.lines.lines) > 0:
+            self.my_logger.debug(f"\n\tLines: {[l.label for l in self.lines.lines]}")
 
     def get_radec_position_after_pm(self, date_obs):
         if self.simbad_table is not None:

spectractor/fit/fit_multispectra.py

@@ -1,5 +1,6 @@
 import matplotlib.pyplot as plt
 from matplotlib import cm
+from matplotlib import colormaps
 from matplotlib import colors
 from matplotlib.ticker import MaxNLocator
 from astropy.io import ascii
@@ -38,7 +39,9 @@ def _build_sim_sample(spectra, aerosols=0.05, ozone=300, pwv=5, angstrom_exponen
 
     Examples
     --------
-    >>> _build_sim_sample([Spectrum("./tests/data/reduc_20170530_134_spectrum.fits")])
+    >>> spectra = _build_sim_sample([Spectrum("./tests/data/reduc_20170530_134_spectrum.fits")])
+    >>> len(spectra)
+    1
     """
     sim_spectra = []
     for spec in spectra:
@@ -656,9 +659,9 @@ def plot_fit(self):  # pragma: no cover
         (array(...
         >>> w.plot_fit()
         """
-        cmap_bwr = copy.copy(cm.get_cmap('bwr'))
+        cmap_bwr = copy.copy(colormaps.get_cmap('bwr'))
         cmap_bwr.set_bad(color='lightgrey')
-        cmap_viridis = copy.copy(cm.get_cmap('viridis'))
+        cmap_viridis = copy.copy(colormaps.get_cmap('viridis'))
         cmap_viridis.set_bad(color='lightgrey')
 
         data = copy.deepcopy(self.data)

spectractor/fit/fit_spectrogram.py

@@ -103,7 +103,10 @@ def __init__(self, spectrum, atmgrid_file_name="", fit_angstrom_exponent=False,
         self.fixed_psf_params = np.array([0, 1, 2, 3, 4, 5, 6, 9])
         self.atm_params_indices = np.array([params.get_index(label) for label in ["VAOD", "angstrom_exp_log10", "ozone [db]", "PWV [mm]"]])
         # A2 is free only if spectrogram is a simulation or if the order 2/1 ratio is not known and flat
-        params.fixed[params.get_index(f"A{self.diffraction_orders[1]}")] = "A2_T" not in self.spectrum.header  # not self.spectrum.disperser.flat_ratio_order_2over1
+        if "A2" in params.labels:
+            params.fixed[params.get_index(f"A{self.diffraction_orders[1]}")] = "A2_T" not in self.spectrum.header
+        if "A3" in params.labels:
+            params.fixed[params.get_index(f"A{self.diffraction_orders[2]}")] = "A3_T" not in self.spectrum.header
         params.fixed[params.get_index(r"shift_x [pix]")] = True  # Delta x
         params.fixed[params.get_index(r"shift_y [pix]")] = True  # Delta y
         params.fixed[params.get_index(r"angle [deg]")] = True  # angle
@@ -207,14 +210,14 @@ def set_mask(self, params=None):
             profile_params[:, 0] = 1
             profile_params[:, 1] = dispersion_law.real + self.simulation.r0.real
             profile_params[:, 2] += dispersion_law.imag # - self.bgd_width
-            psf_cube = self.spectrum.chromatic_psf.build_psf_cube(self.simulation.pixels, profile_params,
-                                                                  fwhmx_clip=3 * parameters.PSF_FWHM_CLIP,
-                                                                  fwhmy_clip=parameters.PSF_FWHM_CLIP, dtype="float32")
-            psf_cube_masked = self.spectrum.chromatic_psf.get_psf_cube_masked(psf_cube, convolve=True)
+            psf_cube_masked = self.spectrum.chromatic_psf.build_psf_cube_masked(self.simulation.pixels, profile_params,
+                                                                                fwhmx_clip=3 * parameters.PSF_FWHM_CLIP,
+                                                                                fwhmy_clip=parameters.PSF_FWHM_CLIP)
+            psf_cube_masked = self.spectrum.chromatic_psf.convolve_psf_cube_masked(psf_cube_masked)
             # make rectangular mask per wavelength
             self.simulation.boundaries[order], self.simulation.psf_cubes_masked[order] = self.spectrum.chromatic_psf.get_boundaries(psf_cube_masked)
             self.simulation.psf_cube_sparse_indices[order], self.simulation.M_sparse_indices[order] = self.spectrum.chromatic_psf.get_sparse_indices(psf_cube_masked)
-        mask = np.sum(self.simulation.psf_cubes_masked[self.diffraction_orders[0]].reshape(psf_cube.shape[0], self.simulation.pixels[0].size), axis=0) == 0
+        mask = np.sum(self.simulation.psf_cubes_masked[self.diffraction_orders[0]].reshape(psf_cube_masked.shape[0], self.simulation.pixels[0].size), axis=0) == 0
         self.W = np.copy(self.W_before_mask)
         self.W[mask] = 0
         self.mask = list(np.where(mask)[0])