Merge pull request #160 from LSSTDESC/debug

Add Gaia spectra and SED errors + small debugs

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.1.0'
+__version__ = '3.1.1'
 __version_info__ = tuple(map(int, __version__.split('.')))

spectractor/extractor/extractor.py

@@ -237,6 +237,8 @@ def __init__(self, spectrum, amplitude_priors_method="noprior", verbose=False, p
         self.amplitude_priors = np.copy(self.spectrum.data)
         if self.amplitude_priors_method == "spectrum":
             self.amplitude_priors_cov_matrix = np.copy(self.spectrum.cov_matrix)
+        else:
+            self.amplitude_priors_cov_matrix = np.diag(self.spectrum.err**2)
         if self.spectrum.data.size != self.Nx:  # must rebin the priors
             old_x = np.linspace(0, 1, self.spectrum.data.size)
             new_x = np.linspace(0, 1, self.Nx)
@@ -513,7 +515,7 @@ def simulate(self, *params):
                     #     L = np.linalg.inv(np.linalg.cholesky(M_dot_W_dot_M))
                     #     cov_matrix = L.T @ L
                     # except np.linalg.LinAlgError:
-                    cov_matrix = np.linalg.inv(M_dot_W_dot_M)
+                    cov_matrix = np.linalg.inv(M_dot_W_dot_M.toarray())
                     amplitude_params = cov_matrix @ (M.T @ W_dot_data)
                     if self.amplitude_priors_method == "positive":
                         amplitude_params[amplitude_params < 0] = 0
@@ -952,7 +954,7 @@ def run_ffm_minimisation(w, method="newton", niter=2):
             w.spectrum.header['PIXSHIFT'] = w.params[r"shift_x [pix]"]
             w.spectrum.header['D2CCD'] = w.params[r"D_CCD [mm]"]
             if len(w.diffraction_orders) >= 2:
-                w.spectrum.header['A2_FIT'] = w.params.values[w.diffraction_orders[1]]
+                w.spectrum.header['A2_FIT'] = w.params[f"A{w.diffraction_orders[1]}"]
             w.spectrum.header["ROTANGLE"] = w.params[r"angle [deg]"]
             w.spectrum.header["AM_FIT"] = w.params["z"]
             # Compute next order contamination

spectractor/fit/fit_spectrogram.py

@@ -15,6 +15,10 @@
 from spectractor.simulation.atmosphere import Atmosphere, AtmosphereGrid
 from spectractor.fit.fitter import (FitWorkspace, FitParameters, run_minimisation, run_minimisation_sigma_clipping,
                                     write_fitparameter_json)
+try:
+    from gaiaspec import getGaia
+except ModuleNotFoundError:
+    getGaia = None
 
 plot_counter = 0
 
@@ -57,7 +61,16 @@ def __init__(self, spectrum, atmgrid_file_name="", fit_angstrom_exponent=False,
 
         """
         if not getCalspec.is_calspec(spectrum.target.label):
-            raise ValueError(f"{spectrum.target.label=} must be a CALSPEC star according to getCalspec package.")
+            if getGaia is not None:
+                is_gaiaspec = getGaia.is_gaiaspec(spectrum.target.label)
+                is_gaia_full = False
+                if is_gaiaspec == False:
+                    is_gaia_full = getGaia.is_gaia_full(spectrum.target.label)
+                if not is_gaiaspec:
+                    if not is_gaia_full:
+                        raise ValueError(f"{spectrum.target.label=} must be a CALSPEC or GAIA star.")
+            else:
+                raise ValueError(f"{spectrum.target.label=} must be a CALSPEC star according to getCalspec package.")
         self.spectrum = spectrum
         self.filename = spectrum.filename.replace("spectrum", "spectrogram")
         self.diffraction_orders = np.arange(spectrum.order, spectrum.order + 3 * np.sign(spectrum.order), np.sign(spectrum.order))
@@ -119,7 +132,7 @@ def __init__(self, spectrum, atmgrid_file_name="", fit_angstrom_exponent=False,
         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
         if "A2" in params.labels:
-            params.fixed[params.get_index(f"A{self.diffraction_orders[1]}")] = False #"A2_T" not in self.spectrum.header
+            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:
             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

spectractor/fit/fit_spectrum.py

@@ -14,7 +14,11 @@
                                     write_fitparameter_json)
 from spectractor.tools import plot_spectrum_simple
 
-
+try:
+    from gaiaspec import getGaia
+except ModuleNotFoundError:
+    getGaia = None
+
 class SpectrumFitWorkspace(FitWorkspace):
 
     def __init__(self, spectrum, atmgrid_file_name="", fit_angstrom_exponent=False,
@@ -58,7 +62,16 @@ def __init__(self, spectrum, atmgrid_file_name="", fit_angstrom_exponent=False,
         """
         self.my_logger = set_logger(self.__class__.__name__)
         if not getCalspec.is_calspec(spectrum.target.label):
-            raise ValueError(f"{spectrum.target.label=} must be a CALSPEC star according to getCalspec package.")
+            if getGaia is not None:
+                is_gaiaspec = getGaia.is_gaiaspec(spectrum.target.label)
+                is_gaia_full = False
+                if is_gaiaspec == False:
+                    is_gaia_full = getGaia.is_gaia_full(spectrum.target.label)
+                if not is_gaiaspec:
+                    if not is_gaia_full:
+                        raise ValueError(f"{spectrum.target.label=} must be a CALSPEC or GAIA star.")
+            else:
+                raise ValueError(f"{spectrum.target.label=} must be a CALSPEC star according to getCalspec package.")
         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
@@ -72,7 +85,7 @@ def __init__(self, spectrum, atmgrid_file_name="", fit_angstrom_exponent=False,
         if not fit_angstrom_exponent:
             fixed[3] = True  # angstrom_exponent
         bounds = [(0, 2), (0, 2/parameters.GRATING_ORDER_2OVER1), (0, 1), (0, 3), (100, 700), (0, 20),
-                       (0.1, 10),(p[7] - 5 * parameters.DISTANCE2CCD_ERR, p[7] + 5 * parameters.DISTANCE2CCD_ERR),
+                       (0.1, 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]",
                                           "reso [nm]", r"D_CCD [mm]", r"alpha_pix [pix]", "B"],

spectractor/fit/fitter.py

@@ -1603,7 +1603,7 @@ def run_minimisation(fit_workspace, method="newton", epsilon=None, xtol=1e-4, ft
         x_scale = np.abs(guess)
         x_scale[x_scale == 0] = 0.1
         p = optimize.least_squares(fit_workspace.weighted_residuals, guess, verbose=2, ftol=1e-6, x_scale=x_scale,
-                                   diff_step=0.001, bounds=bounds.T)
+                                   diff_step=0.001, bounds=list(np.array(bounds).T))
         fit_workspace.params.values = p.x  # m.np_values()
         if verbose:
             my_logger.debug(f"\n\t{p}")

spectractor/simulation/simulator.py

@@ -92,11 +92,11 @@ def simulate_without_atmosphere(self, lambdas):
         self.data *= self.throughput.transmission(lambdas)
         self.data *= self.target.sed(lambdas)
         self.err = np.zeros_like(self.data)
-        idx = np.where(self.throughput.transmission(lambdas) > 0)[0]
-        self.err[idx] = self.throughput.transmission_err(lambdas)[idx] / self.throughput.transmission(lambdas)[idx]
-        self.err[idx] *= self.data[idx]
-        idx = np.where(self.throughput.transmission(lambdas) <= 0)[0]
-        self.err[idx] = 1e6 * np.max(self.err)
+        idx = (self.throughput.transmission(lambdas) > 0) & (self.target.sed(lambdas) > 0)
+        self.err[idx] = np.sqrt((self.throughput.transmission_err(lambdas[idx]) / self.throughput.transmission(lambdas[idx]))**2 + (self.target.sed_err(lambdas[idx])/self.target.sed(lambdas[idx]))**2)
+        self.err[idx] *= np.abs(self.data[idx])
+        idx = (self.throughput.transmission(lambdas) <= 0) | (self.target.sed(lambdas) <= 0)
+        self.err[idx] = 10 * np.max(self.err)
         return self.data, self.err
 
     def simulate(self, A1=1.0, A2=0., aerosols=0.05, angstrom_exponent=None, ozone=300, pwv=5, reso=0.,
@@ -383,10 +383,16 @@ def integrand(lbda):
         spectrum_err = np.zeros_like(spectrum)
         idx = telescope_transmission > 0
         spectrum_err[idx] = self.throughput.transmission_err(lambdas)[idx] / telescope_transmission[idx] * spectrum[idx]
-        # idx = telescope_transmission <= 0: not ready yet to be implemented
+        idx = (telescope_transmission > 0) & (self.target.sed(lambdas) > 0)
+        spectrum_err[idx] = np.sqrt((self.throughput.transmission_err(lambdas[idx]) / telescope_transmission[idx])**2 + (self.target.sed_err(lambdas[idx])/self.target.sed(lambdas[idx]))**2)
+        spectrum_err[idx] *= np.abs(spectrum[idx])
+        idx = (telescope_transmission <= 0) | (self.target.sed(lambdas) <= 0)
+        spectrum_err[idx] = 10 * np.max(spectrum_err)
+        ####TOCHECK idx = telescope_transmission <= 0: not ready yet to be implemented
         # spectrum_err[idx] = 1e6 * np.max(spectrum_err)
         return spectrum, spectrum_err
 
+
     def simulate(self, A1=1.0, A2=0., A3=0., aerosols=0.05, angstrom_exponent=None, ozone=300, pwv=5,
                  D=parameters.DISTANCE2CCD, shift_x=0., shift_y=0., angle=0., psf_poly_params=None):
         """