Test ccdgains

config

@@ -0,0 +1 @@
+./spectractor/config

pyproject.toml

@@ -38,12 +38,18 @@ dependencies = [
 Homepage = "https://github.com/LSSTDESC/Spectractor"
 Repository = "https://github.com/LSSTDESC/Spectractor"
 
+[tool.setuptools]
+packages = [
+    "spectractor",
+    "spectractor.extractor",
+    "spectractor.fit",
+    "spectractor.simulation",
+]
+include-package-data = true
+
 [tool.setuptools.dynamic]
 version = {attr = "spectractor._version.__version__"}
 
-[tool.setuptools.packages.find]
-include = ["spectractor*"]
-
 [tool.setuptools.package-data]
 spectractor = [
     "config/*.ini",
@@ -71,17 +77,12 @@ docs = [
     "sphinxcontrib-apidoc"
 ]
 
-[tool.setuptools.exclude-package-data]
-"*" = ["tests/*"]
-
 [tool.pytest.ini_options]
 testpaths = ["tests"]
 python_files = ["test_*.py"]
 python_classes = ["Test*"]
 python_functions = ["test_*"]
-# Exclure parameters.py et mcmc.py de la collecte de tests
-addopts = "--ignore=spectractor/parameters.py --ignore=spectractor/fit/mcmc.py --ignore=spectractor/fit/astrometry.py"
-norecursedirs = ["spectractor/fit"]
+addopts = "--ignore=spectractor/parameters.py --ignore=spectractor/fit/mcmc.py"
 
 [tool.coverage.run]
 omit = [
@@ -102,4 +103,5 @@ exclude_lines = [
     "if 0:",
     "if False:",
     "if TYPE_CHECKING:",
-]
+]
+show_missing = true

spectractor/extractor/targets.py

@@ -19,14 +19,19 @@
 from spectractor.config import set_logger
 from spectractor.extractor.spectroscopy import (Lines, HGAR_LINES, HYDROGEN_LINES, ATMOSPHERIC_LINES,
                                                 ISM_LINES, STELLAR_LINES)
+
 from getCalspec import getCalspec
 
 # Astroquery versions change the Simbad API.
 _astroquery_version = packaging.version.parse(importlib.metadata.version("astroquery"))
 if _astroquery_version < packaging.version.parse("0.4.8"):
     _USE_NEW_SIMBAD = False
+    _SIMBAD_VOTABLE_FIELDS = ('flux(U)', 'flux(B)', 'flux(V)', 'flux(R)', 'flux(I)', 'flux(J)', 'sptype',
+                              'parallax', 'pm', 'z_value', "IDS", "ids")
 else:
     _USE_NEW_SIMBAD = True
+    _SIMBAD_VOTABLE_FIELDS = ('U', 'B', 'V', 'R', 'I', 'J', 'sp_type',
+                              'parallax', 'propermotions', 'rvz_redshift', "IDS", "ids")                
 
 try:
     from gaiaspec import getGaia
@@ -310,14 +315,7 @@ def load(self):
             simbadQuerier = SimbadClass()
             patchSimbadURL(simbadQuerier)
 
-            if _USE_NEW_SIMBAD:
-                simbadQuerier.add_votable_fields('U', 'B', 'V', 'R', 'I', 'J', 'sp_type',
-                                                 'parallax', 'propermotions', 'rvz_redshift')
-            else:
-                simbadQuerier.add_votable_fields(
-                    'flux(U)', 'flux(B)', 'flux(V)', 'flux(R)', 'flux(I)', 'flux(J)', 'sptype',
-                    'parallax', 'pm', 'z_value'
-                )
+            simbadQuerier.add_votable_fields(*_SIMBAD_VOTABLE_FIELDS)
             self.my_logger.debug(f"\n\tDownload {self.label} coordinates from Simbad...")
             self.simbad_table = simbadQuerier.query_object(astroquery_label)
             self.simbad_table.write(os.path.join(cache_location,f"{cache_file}.ecsv"), overwrite=True)

spectractor/fit/fit_multispectra.py

@@ -303,7 +303,7 @@ def _prepare_data(self):
                 data = []
                 for i in range(1, lambdas_bin_edges.size):
                     lbdas = np.arange(self.lambdas_bin_edges[i - 1], self.lambdas_bin_edges[i] + 1, 1)
-                    data.append(np.trapz(data_func(lbdas), x=lbdas) / self.bin_widths)
+                    data.append(np.trapezoid(data_func(lbdas), x=lbdas) / self.bin_widths)
                 self.data[k] = np.copy(data)
         else:
             for k in range(self.nspectra):
@@ -317,7 +317,7 @@ def _prepare_data(self):
                 data = []
                 for i in range(1, lambdas_bin_edges.size):
                     lbdas = np.arange(self.lambdas_bin_edges[i - 1], self.lambdas_bin_edges[i] + 1, 1)
-                    data.append(np.trapz(data_func(lbdas), x=lbdas) / self.bin_widths)
+                    data.append(np.trapezoid(data_func(lbdas), x=lbdas) / self.bin_widths)
                 self.ref_spectrum_cube.append(np.copy(data))
         else:
             for k in range(self.nspectra):
@@ -337,7 +337,7 @@ def _prepare_data(self):
                         err.append(np.nan)
                     else:
                         lbdas = np.arange(self.lambdas_bin_edges[i - 1], self.lambdas_bin_edges[i] + 1, 1)
-                        err.append(np.sqrt(np.abs(np.trapz(err_func(lbdas), x=lbdas) / self.bin_widths)))
+                        err.append(np.sqrt(np.abs(np.trapezoid(err_func(lbdas), x=lbdas) / self.bin_widths)))
                 self.err[k] = np.copy(err)
         else:
             for k in range(self.nspectra):
@@ -445,7 +445,7 @@ def get_truth(self):
             if self.bin_widths > 0:
                 for i in range(1, self.lambdas_bin_edges.size):
                     lbdas = np.arange(self.lambdas_bin_edges[i - 1], self.lambdas_bin_edges[i] + 1, 1)
-                    self.true_atmospheric_transmission.append(np.trapz(tatm(lbdas), x=lbdas) / self.bin_widths)
+                    self.true_atmospheric_transmission.append(np.trapezoid(tatm(lbdas), x=lbdas) / self.bin_widths)
                     # self.true_atmospheric_transmission.append(quad(tatm, self.lambdas_bin_edges[i - 1],
                     #                                                self.lambdas_bin_edges[i])[0] / self.bin_widths)
             else:
@@ -459,7 +459,7 @@ def get_truth(self):
         if self.bin_widths > 0:
             for i in range(1, self.lambdas_bin_edges.size):
                 lbdas = np.arange(self.lambdas_bin_edges[i - 1], self.lambdas_bin_edges[i] + 1, 1)
-                self.true_instrumental_transmission.append(np.trapz(tinst(lbdas), x=lbdas) / self.bin_widths)
+                self.true_instrumental_transmission.append(np.trapezoid(tinst(lbdas), x=lbdas) / self.bin_widths)
                 # self.true_instrumental_transmission.append(quad(tinst, self.lambdas_bin_edges[i - 1],
                 #                                                 self.lambdas_bin_edges[i])[0] / self.bin_widths)
         else:
@@ -537,7 +537,7 @@ def simulate(self, aerosols, angstrom_exponent, ozone, pwv, reso, *A1s):
                         delta = self.lambdas_bin_edges[i] - self.lambdas_bin_edges[i-1]
                         if delta > 0:
                             lbdas = np.arange(self.lambdas_bin_edges[i-1] + deltas[k], self.lambdas_bin_edges[i] + deltas[k] + 1, 1)
-                            atm.append(np.trapz(a(lbdas), x=lbdas)/delta)
+                            atm.append(np.trapezoid(a(lbdas), x=lbdas)/delta)
                         else:
                             atm.append(1)
                 else:
@@ -717,7 +717,7 @@ def plot_transmissions(self):  # pragma: no cover
         tatm_binned = []
         for i in range(1, self.lambdas_bin_edges.size):
             lbdas = np.arange(self.lambdas_bin_edges[i - 1], self.lambdas_bin_edges[i] + 1, 1)
-            tatm_binned.append(np.trapz(tatm(lbdas), x=lbdas) / (self.lambdas_bin_edges[i] - self.lambdas_bin_edges[i - 1]))
+            tatm_binned.append(np.trapezoid(tatm(lbdas), x=lbdas) / (self.lambdas_bin_edges[i] - self.lambdas_bin_edges[i - 1]))
 
         ax[0, 1].errorbar(self.lambdas[0], tatm_binned,
                           label=r'$T_{\mathrm{atm}}$', fmt='k.')  # , markersize=0.1)
@@ -798,7 +798,7 @@ def save_transmissions(self):
         tatm_binned = []
         for i in range(1, self.lambdas_bin_edges.size):
             lbdas = np.arange(self.lambdas_bin_edges[i - 1], self.lambdas_bin_edges[i] + 1, 1)
-            tatm_binned.append(np.trapz(tatm(lbdas), x=lbdas) / (self.lambdas_bin_edges[i] - self.lambdas_bin_edges[i - 1]))
+            tatm_binned.append(np.trapezoid(tatm(lbdas), x=lbdas) / (self.lambdas_bin_edges[i] - self.lambdas_bin_edges[i - 1]))
         throughput = self.amplitude_params / self.spectra[0].disperser.transmission(self.lambdas[0])
         throughput_err = self.amplitude_params_err / self.spectra[0].disperser.transmission(self.lambdas[0])
         file_name = self.output_file_name + f"_transmissions.txt"

spectractor/fit/fit_spectrogram.py

@@ -103,7 +103,7 @@ def __init__(self, spectrum, atmgrid_file_name="", fit_angstrom_exponent=False,
                        r"$P_{\mathrm{atm}}$ [hPa]"]
         for order in self.diffraction_orders:
             axis_names += [label+rf"$\!_{order}$" for label in psf_poly_params_names]
-        bounds = [[0, 2], [0, 2], [0, 2], [0, 10], [0, 3], [100, 700], [0, 20], [0.8, 1.2], [0, np.inf],
+        bounds = [[0, 2], [0, 2], [0, 2], [0, 10], [0, 4], [100, 700], [0, 20], [0.8, 1.2], [0, np.inf],
                   [D2CCD - 5 * parameters.DISTANCE2CCD_ERR, D2CCD + 5 * parameters.DISTANCE2CCD_ERR], [-2, 2],
                   [-10, 10], [-90, 90], [0, np.inf]]
         bounds += list(psf_poly_params_bounds) * len(self.diffraction_orders)
@@ -132,6 +132,7 @@ def __init__(self, spectrum, atmgrid_file_name="", fit_angstrom_exponent=False,
         params.fixed[params.get_index(f"A{self.diffraction_orders[0]}")] = False  # A1
         self.atm_params_indices = np.array([params.get_index(label) for label in ["VAOD", "angstrom_exp", "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[0]}")] = False  # A1
         if "A2" in params.labels:
             params.fixed[params.get_index(f"A{self.diffraction_orders[1]}")] = False  #not getCalspec.is_calspec(spectrum.target.label) #"A2_T" not in self.spectrum.header
         if "A3" in params.labels:
@@ -609,7 +610,7 @@ def run_spectrogram_minimisation(fit_workspace, method="newton", verbose=False):
         fit_workspace.spectrogram_simulation.fast_sim = False
         fit_workspace.spectrogram_simulation.fix_psf_cube = False
         fit_workspace.params.fixed = [True] * len(fit_workspace.params.values)
-        # fit_workspace.params.fixed[fit_workspace.params.get_index(r"A1")] = False  # A1
+        fit_workspace.params.fixed[fit_workspace.params.get_index(r"A1")] = False  # A1
         fit_workspace.params.fixed[fit_workspace.params.get_index(r"shift_y [pix]")] = False  # shift y
         fit_workspace.params.fixed[fit_workspace.params.get_index(r"angle [deg]")] = False  # angle
         run_minimisation(fit_workspace, "newton", xtol=1e-2, ftol=0.01, with_line_search=False)
@@ -618,19 +619,23 @@ def run_spectrogram_minimisation(fit_workspace, method="newton", verbose=False):
         fit_workspace.spectrogram_simulation.fast_sim = False
         fit_workspace.spectrogram_simulation.fix_psf_cube = False
         fit_workspace.params.fixed = np.copy(fixed_default)
+        # fit_workspace.params.values[fit_workspace.params.get_index(r"A1")] = 1.0  # A1
         # guess = fit_workspace.p
         # params_table, costs = run_gradient_descent(fit_workspace, guess, epsilon, params_table, costs,
         #                                            fix=fit_workspace.fixed, xtol=1e-6, ftol=1 / fit_workspace.data.size,
         #                                            niter=40)
         run_minimisation_sigma_clipping(fit_workspace, method="newton", xtol=1e-6,
                                         ftol=1 / fit_workspace.data.size, sigma_clip=100, niter_clip=3, verbose=verbose,
                                         with_line_search=True)
+        extra = {"chi2": fit_workspace.costs[-1] / fit_workspace.data.size,
+                 "date-obs": fit_workspace.spectrum.date_obs,
+                 "outliers": len(fit_workspace.outliers)}
+        fit_workspace.params.extra = extra
         my_logger.info(f"\n\tNewton: total computation time: {time.time() - start}s")
         if fit_workspace.filename != "":
             fit_workspace.params.plot_correlation_matrix()
             write_fitparameter_json(fit_workspace.params.json_filename, fit_workspace.params,
-                                    extra={"chi2": fit_workspace.costs[-1] / fit_workspace.data.size,
-                                           "date-obs": fit_workspace.spectrum.date_obs})
+                                    extra=extra)
             # save_gradient_descent(fit_workspace, costs, params_table)
             fit_workspace.plot_fit()
 

spectractor/fit/fit_spectrum.py

@@ -75,16 +75,16 @@ def __init__(self, spectrum, atmgrid_file_name="", fit_angstrom_exponent=False,
         self.spectrum = spectrum
         p = np.array([1, 0, 0.05, 1.2, 400, 5, 1, self.spectrum.header['D2CCD'], self.spectrum.header['PIXSHIFT'], 0])
         fixed = [False] * p.size
-        fixed[0] = True
-        fixed[1] = "A2_T" not in self.spectrum.header  # fit A2 only on sims to evaluate extraction biases
+        fixed[0] = False
+        fixed[1] = False #"A2_T" not in self.spectrum.header  # fit A2 only on sims to evaluate extraction biases
         fixed[5] = False
         # fixed[6:8] = [True, True]
         fixed[8] = True
         fixed[9] = True
         # fixed[-1] = True
         if not fit_angstrom_exponent:
             fixed[3] = True  # angstrom_exponent
-        bounds = [(0, 2), (0, 2/parameters.GRATING_ORDER_2OVER1), (0, 10), (0, 3), (100, 700), (0, 20),
+        bounds = [(0, 2), (0, 2/parameters.GRATING_ORDER_2OVER1), (0, 10), (0, 4), (100, 700), (0, 20),
                        (0.5, 20),(p[7] - 5 * parameters.DISTANCE2CCD_ERR, p[7] + 5 * parameters.DISTANCE2CCD_ERR),
                   (-2, 2), (-np.inf, np.inf)]
         params = FitParameters(p, labels=["A1", "A2", "VAOD", "angstrom_exp", "ozone [db]", "PWV [mm]",
@@ -396,7 +396,8 @@ def run_spectrum_minimisation(fit_workspace, method="newton", sigma_clip=20):
 
         fit_workspace.simulation.fast_sim = False
         fixed = copy.copy(fit_workspace.params.fixed)
-        fit_workspace.params.fixed[6] = True
+        fit_workspace.params.fixed = [True] * len(fit_workspace.params)
+        fit_workspace.params.fixed[0] = False
         run_minimisation(fit_workspace, method="newton", xtol=1e-3, ftol=100 / fit_workspace.data.size,
                          verbose=False)
         fit_workspace.params.fixed = fixed
@@ -405,10 +406,12 @@ def run_spectrum_minimisation(fit_workspace, method="newton", sigma_clip=20):
 
         fit_workspace.params.plot_correlation_matrix()
         fit_workspace.plot_fit()
+        extra = {"chi2": fit_workspace.costs[-1] / fit_workspace.data.size,
+                 "date-obs": fit_workspace.spectrum.date_obs,
+                 "outliers": len(fit_workspace.outliers)}
+        fit_workspace.params.extra = extra
         if fit_workspace.filename != "":
-            write_fitparameter_json(fit_workspace.params.json_filename, fit_workspace.params,
-                                    extra={"chi2": fit_workspace.costs[-1] / fit_workspace.data.size,
-                                           "date-obs": fit_workspace.spectrum.date_obs})
+            write_fitparameter_json(fit_workspace.params.json_filename, fit_workspace.params, extra=extra)
             # save_gradient_descent(fit_workspace, costs, params_table)
 
 

spectractor/fit/fitter.py

@@ -1328,12 +1328,12 @@ def gradient_descent(fit_workspace, niter=10, xtol=1e-3, ftol=1e-3, with_line_se
             if ip not in ipar:
                 continue
             # check for null vectors
-            if J_norms[ip] < 1e-20:
+            if J_norms[ip] < np.sqrt(len(J_vectors[ip]))*np.finfo(np.float64).eps and len(np.where(J_vectors[ip]==0)[0]) > J_vectors[ip].size // 2:
                 ipar = np.delete(ipar, list(ipar).index(ip))
                 fit_workspace.params.fixed[ip] = True
                 my_logger.warning(
                     f"\n\tStep {i}: {fit_workspace.params.labels[ip]} has a null Jacobian; parameter is fixed "
-                    f"at its last known current value ({tmp_params[ip]}).")
+                    f"at its last known current value ({tmp_params[ip]}) because |J[par]|={J_norms[ip]:.3e}<{np.sqrt(J.shape[1])*np.finfo(np.float64).eps} with more than half values being zeros.")
                 continue
             # check for degeneracies using Cauchy-Schwartz inequality; fix the second parameter
             for jp in range(ip, J.shape[0]):

spectractor/parameters.py

@@ -48,7 +48,7 @@ def __getattr__(name):
 
 # Paths
 DISPERSER_DIR = "./extractor/dispersers/"
-CONFIG_DIR = "../config/"
+CONFIG_DIR = "./config/"
 THROUGHPUT_DIR = "./simulation/CTIOThroughput/"
 if 'ASTROMETRYNET_DIR' in os.environ:
     ASTROMETRYNET_DIR = os.getenv('ASTROMETRYNET_DIR') + '/'