Merge pull request #178 from LSSTDESC/test_ccdgains_robustfit

Test ccdgains robustfit

Author: jeremyneveu

pyproject.toml

@@ -39,24 +39,25 @@ 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.packages.find]
+where = ["."]
+include = ["spectractor*"]
+
 [tool.setuptools.dynamic]
 version = {attr = "spectractor._version.__version__"}
 
 [tool.setuptools.package-data]
+"*" = ["*.ini", "*.txt"]
 spectractor = [
     "config/*.ini",
-    "extractor/dispersers/**/*",
-    "simulation/CTIOThroughput/*.txt",
-    "simulation/AuxTelThroughput/*.txt", 
-    "simulation/StarDiceThroughput/*.txt"
+    "extractor/dispersers/**/*"
+]
+"spectractor.simulation" = [
+    "CTIOThroughput/*.txt",
+    "AuxTelThroughput/*.txt",
+    "StarDiceThroughput/*.txt"
 ]
 
 [project.optional-dependencies]

spectractor/extractor/extractor.py

@@ -495,7 +495,7 @@ def simulate(self, *params):
         if self.flat is not None:
             # multiply each M matrix columns by the flat array (see the docstring)
             # TODO: if flat array is a cube flat, needs to multiply directly in build_sparse_M
-            dia = sparse.dia_matrix(([self.flat], [0]), shape=(self.flat.size, self.flat.size))
+            dia = sparse.dia_matrix(([self.flat.astype("float32")], [0]), shape=(self.flat.size, self.flat.size))
             M = (dia @ M).tocsc()
 
 

spectractor/fit/fit_spectrogram.py

@@ -104,7 +104,7 @@ def __init__(self, spectrum, atmgrid_file_name="", fit_angstrom_exponent=False,
         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, 4], [100, 700], [0, 20], [0.8, 1.2], [0, np.inf],
-                  [D2CCD - 5 * parameters.DISTANCE2CCD_ERR, D2CCD + 5 * parameters.DISTANCE2CCD_ERR], [-2, 2],
+                  [D2CCD - 20 * parameters.DISTANCE2CCD_ERR, D2CCD + 20 * parameters.DISTANCE2CCD_ERR], [-10, 10],
                   [-10, 10], [-90, 90], [0, np.inf]]
         bounds += list(psf_poly_params_bounds) * len(self.diffraction_orders)
         fixed = [False] * p.size
@@ -132,9 +132,9 @@ 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
+        params.fixed[params.get_index(f"A{self.diffraction_orders[0]}")] = True  # 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
+            params.fixed[params.get_index(f"A{self.diffraction_orders[1]}")] = True  #not getCalspec.is_calspec(spectrum.target.label) #"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]")] = False  # Delta x
@@ -181,6 +181,11 @@ def __init__(self, spectrum, atmgrid_file_name="", fit_angstrom_exponent=False,
         if not fit_angstrom_exponent:
             self.params.fixed[self.params.get_index("angstrom_exp")] = True  # angstrom exponent
         self.params.values[self.params.get_index("angstrom_exp")] = self.atmosphere.angstrom_exponent_default
+        if np.min(self.spectrum.lambdas) > 500:
+            self.fit_angstrom_exponent = False
+            self.params.fixed[self.params.get_index("angstrom_exp")] = True
+            self.params.values[self.params.get_index("angstrom_exp")] = 0
+            self.my_logger.warning("\n\tWavelengths below 500nm detected: angstrom exponent fitting disabled and fixed to 0.")
         self.spectrogram_simulation = SpectrogramModel(self.spectrum, atmosphere=self.atmosphere,
                                            diffraction_orders=self.diffraction_orders,
                                            fast_sim=False, with_adr=True)
@@ -577,7 +582,7 @@ def run_spectrogram_minimisation(fit_workspace, method="newton", verbose=False):
     my_logger = set_logger(__name__)
     guess = np.asarray(fit_workspace.params.values)
     fit_workspace.simulate(*guess)
-    fit_workspace.plot_fit()
+    # fit_workspace.plot_fit()
     if method != "newton":
         run_minimisation(fit_workspace, method=method)
     else:
@@ -610,7 +615,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"VAOD")] = False  # VAOD
         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)
@@ -624,8 +629,8 @@ def run_spectrogram_minimisation(fit_workspace, method="newton", verbose=False):
         # 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,
+        run_minimisation_sigma_clipping(fit_workspace, method="newton", xtol=1e-10,
+                                        ftol=1e-3 / 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,

spectractor/fit/fit_spectrum.py

@@ -75,18 +75,18 @@ 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] = False
-        fixed[1] = False #"A2_T" not in self.spectrum.header  # fit A2 only on sims to evaluate extraction biases
+        fixed[0] = True  # A1
+        fixed[1] = True  #"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[8] = False  # alpha_pix
         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, 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)]
+                  (0.5, 20),(p[7] - 20 * parameters.DISTANCE2CCD_ERR, p[7] + 20 * parameters.DISTANCE2CCD_ERR),
+                  (-10, 10), (-np.inf, np.inf)]
         params = FitParameters(p, labels=["A1", "A2", "VAOD", "angstrom_exp", "ozone [db]", "PWV [mm]",
                                           "reso [nm]", r"D_CCD [mm]", r"alpha_pix [pix]", "B"],
                                axis_names=["$A_1$", "$A_2$", "VAOD", r'$\"a$', "ozone [db]", "PWV [mm]",
@@ -112,10 +112,15 @@ def __init__(self, spectrum, atmgrid_file_name="", fit_angstrom_exponent=False,
             self.params.fixed[self.params.get_index("angstrom_exp")] = True  # angstrom exponent
             if parameters.VERBOSE:
                 self.my_logger.info(f'\n\tUse atmospheric grid models from file {atmgrid_file_name}. ')
-        self.params.values[self.params.get_index("angstrom_exp")] = self.atmosphere.angstrom_exponent_default
+
         self.lambdas = self.spectrum.lambdas
         self.fit_angstrom_exponent = fit_angstrom_exponent
         self.params.values[self.params.get_index("angstrom_exp")] = self.atmosphere.angstrom_exponent_default
+        if np.min(self.spectrum.lambdas) > 500:
+            self.fit_angstrom_exponent = False
+            self.params.fixed[self.params.get_index("angstrom_exp")] = True
+            self.params.values[self.params.get_index("angstrom_exp")] = 0
+            self.my_logger.warning("\n\tWavelengths below 500nm detected: angstrom exponent fitting disabled and fixed to 0.")
         self.simulation = SpectrumSimulation(self.spectrum, atmosphere=self.atmosphere, fast_sim=True, with_adr=True)
         self.amplitude_truth = None
         self.lambdas_truth = None
@@ -397,12 +402,27 @@ 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 = [True] * len(fit_workspace.params)
-        fit_workspace.params.fixed[0] = False
+        fit_workspace.params.fixed[fit_workspace.params.get_index(r"VAOD")] = False
         run_minimisation(fit_workspace, method="newton", xtol=1e-3, ftol=100 / fit_workspace.data.size,
-                         verbose=False)
+                         verbose=False, with_line_search=False)
         fit_workspace.params.fixed = fixed
         run_minimisation_sigma_clipping(fit_workspace, method="newton", xtol=1e-6,
-                                        ftol=1 / fit_workspace.data.size, sigma_clip=sigma_clip, niter_clip=3, verbose=False)
+                                        ftol=1e-3 / fit_workspace.data.size, sigma_clip=sigma_clip, niter_clip=3, verbose=False)
+
+        # alternate fixing dccd and pixshift with fitting all parameters
+        for i in range(3):
+            fixed = copy.copy(fit_workspace.params.fixed)
+            fit_workspace.params.fixed = [True] * len(fit_workspace.params)
+            fit_workspace.params.fixed[6] = False  # reso
+            fit_workspace.params.fixed[7] = False  # dccd
+            fit_workspace.params.fixed[8] = False  # pixshift
+            run_minimisation_sigma_clipping(fit_workspace, method="newton", xtol=1e-6,
+                                        ftol=1e-3 / fit_workspace.data.size, sigma_clip=sigma_clip, niter_clip=3, verbose=False, with_line_search=True)
+            fit_workspace.params.fixed = fixed
+            run_minimisation_sigma_clipping(fit_workspace, method="newton", xtol=1e-6,
+                                        ftol=1e-3 / fit_workspace.data.size, sigma_clip=sigma_clip, niter_clip=3, verbose=False, with_line_search=True)
+            
+         
 
         fit_workspace.params.plot_correlation_matrix()
         fit_workspace.plot_fit()

spectractor/fit/fitter.py

@@ -1277,6 +1277,91 @@ def save_gradient_descent(self):
         np.savetxt(output_filename, t.T, header=h, delimiter=",")
         self.my_logger.info(f"\n\tSave gradient descent log {output_filename}.")
 
+    def plot_triangle_fisher(self, nsamples=10000, max_params=8):
+        """Plot a triangle (corner) plot of probability contours and distributions using Fisher matrix method.
+        
+        This method generates samples from a multivariate normal distribution based on the parameter 
+        covariance matrix and creates a corner plot showing 1D and 2D marginal distributions.
+        
+        Parameters
+        ----------
+        nsamples: int, optional
+            Number of samples to generate from the multivariate normal distribution (default: 10000).
+        max_params: int, optional
+            Maximum number of parameters to display. Only the first max_params free parameters 
+            will be shown (default: 8).
+            
+        Examples
+        --------
+        >>> from spectractor.fit.fitter import FitParameters, FitWorkspace
+        >>> params = FitParameters(values=[1, 2, 3], axis_names=["x", "y", "z"])
+        >>> params.cov = np.array([[1, -0.5, 0], [-0.5, 1, -1], [0, -1, 1]])
+        >>> w = FitWorkspace(params=params, filename="test.fits")
+        >>> w.plot_triangle_fisher(nsamples=1000, max_params=3)
+        
+        Notes
+        -----
+        This method requires the `corner` package to be installed.
+        The plot is saved if parameters.SAVE is True with suffix '_triangle.pdf'.
+        """
+        try:
+            import corner
+        except ImportError:
+            self.my_logger.warning("Package 'corner' not installed. Cannot create triangle plot. "
+                                   "Install it with: pip install corner")
+            return
+        
+        # Get free parameters indices
+        ipar = self.params.get_free_parameters()
+        
+        # Limit to max_params parameters
+        n_params = min(len(ipar), max_params)
+        ipar = ipar[:n_params]
+        
+        if len(ipar) == 0:
+            self.my_logger.warning("No free parameters to plot.")
+            return
+            
+        # Check that covariance matrix is available
+        if self.params.cov.size == 0:
+            self.my_logger.warning("Covariance matrix not available. Run fit first.")
+            return
+        
+        # Extract mean values and covariance matrix for free parameters
+        mean_values = self.params.values[ipar]
+        cov_matrix = self.params.cov[:n_params, :n_params]
+        
+        # Check if covariance matrix is positive definite
+        eigenvalues = np.linalg.eigvals(cov_matrix)
+        if np.any(eigenvalues <= 0):
+            self.my_logger.warning("Covariance matrix is not positive definite. "
+                                   "Cannot generate samples.")
+            return
+        
+        # Generate samples from multivariate normal distribution
+        samples = np.random.multivariate_normal(mean_values, cov_matrix, size=nsamples)
+        
+        # Get axis names for the selected parameters
+        axis_names = [self.params.axis_names[i] for i in ipar]
+        
+        # Create corner plot
+        fig = corner.corner(samples, labels=axis_names, 
+                           quantiles=[0.16, 0.5, 0.84],
+                           show_titles=True,
+                           title_kwargs={"fontsize": 12},
+                           truths=mean_values)
+        
+        # Save figure if requested
+        if (parameters.SAVE or parameters.LSST_SAVEFIGPATH) and self.filename != "":  # pragma: no cover
+            figname = os.path.splitext(self.filename)[0] + "_triangle.pdf"
+            self.my_logger.info(f"\n\tSave triangle plot {figname}.")
+            fig.savefig(figname, dpi=100, bbox_inches='tight')
+        
+        if parameters.DISPLAY:  # pragma: no cover
+            plt.show()
+        
+        return fig
+
 
 def gradient_descent(fit_workspace, niter=10, xtol=1e-3, ftol=1e-3, with_line_search=True):
     """
@@ -1333,7 +1418,7 @@ def gradient_descent(fit_workspace, niter=10, xtol=1e-3, ftol=1e-3, with_line_se
                 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]}) because |J[par]|={J_norms[ip]:.3e}<{np.sqrt(J.shape[1])*np.finfo(np.float64).eps} with more than half values being zeros.")
+                    f"at its last known current value ({tmp_params[ip]}) because |J[par]|={J_norms[ip]:.3e}<{np.sqrt(len(J_vectors[ip]))*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]):