Merge pull request #140 from LSSTDESC/v3.0.2

V3.0.2

config/auxtel.ini

@@ -15,6 +15,10 @@ SPECTRACTOR_DECONVOLUTION_PSF2D = True
 SPECTRACTOR_DECONVOLUTION_FFM = True
 # value of sigma clip parameter for the spectractor deconvolution process PSF2D and FFM
 SPECTRACTOR_DECONVOLUTION_SIGMA_CLIP = 100
+# maximum time per gradient descent iteration before TimeoutError in seconds
+SPECTRACTOR_FIT_TIMEOUT_PER_ITER = 600
+# maximum time per gradient descent before TimeoutError in seconds
+SPECTRACTOR_FIT_TIMEOUT = 3600
 
 [instrument]
 # instrument name
@@ -23,20 +27,16 @@ OBS_NAME = AUXTEL
 OBS_ALTITUDE = 2.66299616375123
 # observatory latitude
 OBS_LATITUDE = -30.2446389756252
-# diameter of the telescope in meter
-OBS_DIAMETER = 1.20
+# effective surface of the telescope in cm**2 accounting for obscuration
+OBS_SURFACE = 9636
 # observation epoch
 OBS_EPOCH = J2000.0
-# systematics on the instrument transmission
-OBS_TRANSMISSION_SYSTEMATICS = 0.005
 # observed object to choose between STAR, HG-AR, MONOCHROMATOR
 OBS_OBJECT_TYPE = STAR
-# telescope transmission file
-OBS_TELESCOPE_TRANSMISSION = multispectra_holo4_003_HD142331_AuxTel_throughput.txt
 # full instrument transmission file
-OBS_FULL_INSTRUMENT_TRANSMISSON = multispectra_holo4_003_multi_AuxTel_throughput.txt
-# quantum efficiency of the detector file
-OBS_QUANTUM_EFFICIENCY = multispectra_holo4_003_HD142331_AuxTel_throughput.txt
+OBS_FULL_INSTRUMENT_TRANSMISSON = multispectra_holo4_003_HD142331_AuxTel_throughput.txt
+# systematics on the instrument transmission
+OBS_TRANSMISSION_SYSTEMATICS = 0.005
 # Camera (x,y) rotation angle with respect to (north-up, east-left) system
 OBS_CAMERA_ROTATION = 0
 # Camera (x,y) flip signs with respect to (north-up, east-left) system

config/auxtel_quicklook.ini

@@ -21,20 +21,16 @@ OBS_NAME = AUXTEL
 OBS_ALTITUDE = 2.66299616375123
 # observatory latitude
 OBS_LATITUDE = -30.2446389756252
-# diameter of the telescope in meter
-OBS_DIAMETER = 1.20
+# effective surface of the telescope in cm**2 accounting for obscuration
+OBS_SURFACE = 9636
 # observation epoch
 OBS_EPOCH = J2000.0
-# systematics on the instrument transmission
-OBS_TRANSMISSION_SYSTEMATICS = 0.005
 # observed object to choose between STAR, HG-AR, MONOCHROMATOR
 OBS_OBJECT_TYPE = STAR
-# telescope transmission file
-OBS_TELESCOPE_TRANSMISSION = calexp_2020031500162-EMPTY_ronchi90lpmm-det000_auxtel_transmission.txt
 # full instrument transmission file
-OBS_FULL_INSTRUMENT_TRANSMISSON = calexp_2020031500162-EMPTY_ronchi90lpmm-det000_auxtel_transmission.txt
-# quantum efficiency of the detector file
-OBS_QUANTUM_EFFICIENCY = calexp_2020031500162-EMPTY_ronchi90lpmm-det000_auxtel_transmission.txt
+OBS_FULL_INSTRUMENT_TRANSMISSON = multispectra_holo4_003_HD142331_AuxTel_throughput.txt
+# systematics on the instrument transmission
+OBS_TRANSMISSION_SYSTEMATICS = 0.005
 # Camera (x,y) rotation angle with respect to (north-up, east-left) system
 OBS_CAMERA_ROTATION = 0
 # Camera (x,y) flip signs with respect to (north-up, east-left) system

config/ctio.ini

@@ -23,20 +23,16 @@ OBS_NAME = CTIO
 OBS_ALTITUDE = 2.200
 # observatory latitude
 OBS_LATITUDE = -30 10 07.90
-# diameter of the telescope in meter
-OBS_DIAMETER = 0.9
+# effective surface of the telescope in cm**2 accounting for obscuration
+OBS_SURFACE = 6361
 # observation epoch
 OBS_EPOCH = J2000.0
-# systematics on the instrument transmission
-OBS_TRANSMISSION_SYSTEMATICS = 0.0
 # observed object to choose between STAR, HG-AR, MONOCHROMATOR
 OBS_OBJECT_TYPE = STAR
-# telescope transmission file
-OBS_TELESCOPE_TRANSMISSION = ctio_throughput.txt
 # full instrument transmission file
 OBS_FULL_INSTRUMENT_TRANSMISSON = multispectra_Thor300_HD111980_CTIO_throughput_prod7.5.4.txt
-# quantum efficiency of the detector file
-OBS_QUANTUM_EFFICIENCY = qecurve.txt
+# systematics on the instrument transmission
+OBS_TRANSMISSION_SYSTEMATICS = 0.0
 # Camera (x,y) rotation angle with respect to (north-up, east-left) system
 OBS_CAMERA_ROTATION = 180
 # Camera (x,y) flip signs with respect to (north-up, east-left) system

config/default.ini

@@ -15,6 +15,10 @@ SPECTRACTOR_DECONVOLUTION_PSF2D = True
 SPECTRACTOR_DECONVOLUTION_FFM = True
 # value of sigma clip parameter for the spectractor deconvolution process PSF2D and FFM
 SPECTRACTOR_DECONVOLUTION_SIGMA_CLIP = 20
+# maximum time per gradient descent iteration before TimeoutError in seconds
+SPECTRACTOR_FIT_TIMEOUT_PER_ITER = 600
+# maximum time per gradient descent before TimeoutError in seconds
+SPECTRACTOR_FIT_TIMEOUT = 3600
 
 [instrument]
 # instrument name
@@ -23,20 +27,16 @@ OBS_NAME = DEFAULT
 OBS_ALTITUDE = 0
 # observatory latitude
 OBS_LATITUDE = 0
-# diameter of the telescope in meter
-OBS_DIAMETER = 1
+# effective surface of the telescope in cm**2 accounting for obscuration
+OBS_SURFACE = 1
 # observation epoch
 OBS_EPOCH = J2000.0
-# systematics on the instrument transmission
-OBS_TRANSMISSION_SYSTEMATICS = 0
 # observed object to choose between STAR, HG-AR, MONOCHROMATOR
 OBS_OBJECT_TYPE = STAR
-# telescope transmission file
-OBS_TELESCOPE_TRANSMISSION = ctio_throughput.txt
 # full instrument transmission file
 OBS_FULL_INSTRUMENT_TRANSMISSON = ctio_throughput_300517_v2.txt
-# quantum efficiency of the detector file
-OBS_QUANTUM_EFFICIENCY = qecurve.txt
+# systematics on the instrument transmission
+OBS_TRANSMISSION_SYSTEMATICS = 0
 
 [CCD]
 # size of the image in pixel

config/lpnhe.ini

@@ -11,8 +11,8 @@ OBS_NAME = LPNHE
 OBS_ALTITUDE = 0
 # observatory latitude
 OBS_LATITUDE = 0
-# diameter of the telescope in meter
-OBS_DIAMETER = 0.0762
+# effective surface of the telescope in cm**2 accounting for obscuration
+OBS_SURFACE = 45.60
 # observation epoch
 OBS_EPOCH = J2000.0
 # systematics on the instrument transmission

config/stardice.ini

@@ -26,22 +26,16 @@ OBS_NAME = STARDICE
 OBS_ALTITUDE = 0.650
 # observatory latitude
 OBS_LATITUDE = 43.930832
-# diameter of the telescope in meter
-OBS_DIAMETER = 0.4
-#diameter of the telescope secondary mirror in meter
-OBS_SECONDARY_DIAMETER = 0.11
+# effective surface of the telescope in cm**2 accounting for obscuration
+OBS_SURFACE = 1161.6
 # observation epoch
 OBS_EPOCH = J2000.0
-# systematics on the instrument transmission
-OBS_TRANSMISSION_SYSTEMATICS = 0.005
 # observed object to choose between STAR, HG-AR, MONOCHROMATOR
 OBS_OBJECT_TYPE = STAR
-# telescope transmission file
-OBS_TELESCOPE_TRANSMISSION = StarDice_EMPTY_response_75um_pinhole.txt
 # full instrument transmission file
 OBS_FULL_INSTRUMENT_TRANSMISSON = StarDice_EMPTY_response_75um_pinhole.txt
-# quantum efficiency of the detector file
-OBS_QUANTUM_EFFICIENCY = StarDice_EMPTY_response_75um_pinhole.txt
+# systematics on the instrument transmission
+OBS_TRANSMISSION_SYSTEMATICS = 0.005
 # Camera (x,y) rotation angle with respect to (north-up, east-left) system
 OBS_CAMERA_ROTATION = 180
 # Camera (x,y) flip signs with respect to (north-up, east-left) system

requirements.txt

@@ -20,4 +20,5 @@ deprecated
 pyyaml
 nose
 getCalspec>=2.0.0
-lsst_utils
+lsst_utils
+packaging

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

spectractor/config.py

@@ -131,11 +131,8 @@ def update_derived_parameters():
     parameters.CALIB_BGD_NPARAMS = parameters.CALIB_BGD_ORDER + 1
     parameters.LAMBDAS = np.arange(parameters.LAMBDA_MIN, parameters.LAMBDA_MAX, parameters.LAMBDA_STEP)
     parameters.CCD_ARCSEC2RADIANS = np.pi / (180. * 3600.)  # conversion factor from arcsec to radians
-    if not isinstance(parameters.OBS_DIAMETER, astropy.units.quantity.Quantity):
-        parameters.OBS_DIAMETER = parameters.OBS_DIAMETER * units.m  # Diameter of the telescope
-    if not isinstance(parameters.OBS_SECONDARY_DIAMETER, astropy.units.quantity.Quantity):
-        parameters.OBS_SECONDARY_DIAMETER = parameters.OBS_SECONDARY_DIAMETER * units.m  # Diameter of the telescope
-    parameters.OBS_SURFACE = np.pi * parameters.OBS_DIAMETER ** 2 / 4. - np.pi * parameters.OBS_SECONDARY_DIAMETER ** 2 / 4.  # Surface of telescope
+    if not isinstance(parameters.OBS_SURFACE, astropy.units.quantity.Quantity):
+        parameters.OBS_SURFACE *= units.cm ** 2  # Surface of telescope
     # Conversion factor
     # Units of SEDs in flam (erg/s/cm2/nm) :
     parameters.hc = const.h * const.c  # h.c product of fundamental constants c and h

spectractor/extractor/spectrum.py

@@ -849,12 +849,16 @@ def load_spectrum(self, input_file_name, spectrogram_file_name_override=None,
         # check the version of the file
         if "VERSION" in self.header:
             from spectractor._version import __version__
+            from packaging import version
             if self.config != "":
                 raise AttributeError(f"With Spectractor above 2.4 do not provide a config file in Spectrum(config=...)."
                                      f"Now config parameters are loaded from the file header. Got {self.config=}.")
             if self.header["VERSION"] != str(__version__):
+                self.my_logger.debug(f"\n\tSpectrum file spectractor version {self.header['VERSION']} is "
+                                     f"different from current Spectractor software {__version__}.")
+            if version.parse(self.header["VERSION"]) < version.parse("3.0"):
                 self.my_logger.warning(f"\n\tSpectrum file spectractor version {self.header['VERSION']} is "
-                                       f"different from current Spectractor software {__version__}.")
+                                       f"below Spectractor software 3.0. It may be deprecated.")
             self.load_spectrum_latest(input_file_name)
         else:
             self.my_logger.warning("\n\tNo information about Spectractor software version is given in the header. "

spectractor/fit/fit_spectrogram.py

@@ -85,7 +85,7 @@ def __init__(self, spectrum, atmgrid_file_name="", fit_angstrom_exponent=False,
                        r"$\Delta_{\mathrm{x}}$ [pix]", r"$\Delta_{\mathrm{y}}$ [pix]", r"$\theta$ [deg]", "$B$"]
         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, 0.1], [-5, 2], [100, 700], [0, 10],
+        bounds = [[0, 2], [0, 2], [0, 2], [0, 0.1], [-5, 2], [100, 700], [0, 20],
                   [D2CCD - 5 * parameters.DISTANCE2CCD_ERR, D2CCD + 5 * parameters.DISTANCE2CCD_ERR], [-2, 2],
                   [-10, 10], [-90, 90], [0.8, 1.2]]
         bounds += list(psf_poly_params_bounds) * len(self.diffraction_orders)

spectractor/fit/fit_spectrum.py

@@ -71,7 +71,7 @@ def __init__(self, spectrum, atmgrid_file_name="", fit_angstrom_exponent=False,
         # fixed[-1] = True
         if not fit_angstrom_exponent:
             fixed[3] = True  # angstrom_exponent
-        bounds = [(0, 2), (0, 2/parameters.GRATING_ORDER_2OVER1), (0, 0.1), (-5, 2), (100, 700), (0, 10),
+        bounds = [(0, 2), (0, 2/parameters.GRATING_ORDER_2OVER1), (0, 0.1), (-5, 2), (100, 700), (0, 20),
                        (0.1, 10),(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_log10", "ozone [db]", "PWV [mm]",

spectractor/fit/fitter.py

@@ -1248,6 +1248,7 @@ def gradient_descent(fit_workspace, epsilon, niter=10, xtol=1e-3, ftol=1e-3, wit
     params_table = []
     inv_JT_W_J = np.zeros((len(ipar), len(ipar)), dtype=float)
 
+    start_total = time.time()
     for i in range(niter):
         start = time.time()
         cost, tmp_lambdas, tmp_model, tmp_model_err = fit_workspace.chisq(tmp_params, model_output=True)
@@ -1390,6 +1391,12 @@ def line_search(alpha):
             my_logger.info(f"\n\tGradient descent terminated in {i} iterations because the "
                            f"relative change of cost is below ftol={ftol}.")
             break
+        if time.time() - start > parameters.SPECTRACTOR_FIT_TIMEOUT_PER_ITER:
+            raise TimeoutError(f"Gradient descent iteration longer than {parameters.SPECTRACTOR_FIT_TIMEOUT_PER_ITER=}. "
+                               f"Check data quality or increase parameters.SPECTRACTOR_FIT_TIMEOUT_PER_ITER in config file.")
+    if time.time() - start_total > parameters.SPECTRACTOR_FIT_TIMEOUT:
+        raise TimeoutError(f"Gradient descent longer than {parameters.SPECTRACTOR_FIT_TIMEOUT=}. "
+                           f"Check data quality or increase parameters.SPECTRACTOR_FIT_TIMEOUT in config file.")
     return tmp_params, inv_JT_W_J, np.array(costs), np.array(params_table)
 
 

spectractor/parameters.py

@@ -37,6 +37,8 @@ def __getattr__(name):
 SPECTRACTOR_DECONVOLUTION_FFM = True  # deconvolve spectrogram with full forward model: False, True
 SPECTRACTOR_DECONVOLUTION_SIGMA_CLIP = 20  # value of sigma clip parameter for the spectractor deconvolution process PSF2D and FFM
 SPECTRACTOR_BACKGROUND_SUBTRACTION = True #if True the background is estimated and subtracted 
+SPECTRACTOR_FIT_TIMEOUT_PER_ITER = 600  # maximum time per gradient descent iteration before TimeoutError in seconds
+SPECTRACTOR_FIT_TIMEOUT = 3600  # maximum time per gradient descent before TimeoutError in seconds
 
 # Paths
 DISPERSER_DIR = "./extractor/dispersers/"
@@ -63,14 +65,11 @@ def __getattr__(name):
 OBS_NAME = 'CTIO'
 OBS_ALTITUDE = 2.200  # CTIO altitude in k meters from astropy package (Cerro Pachon)
 OBS_LATITUDE = '-30 10 07.90'  # CTIO latitude
-OBS_DIAMETER = 0.9  # Diameter of the telescope in meter
-OBS_SECONDARY_DIAMETER = 0  # Diameter of the secondary mirror in meter
+OBS_SURFACE = 6361  # Effective surface of the telescope in cm**2 accounting for obscuration
 OBS_EPOCH = "J2000.0"
-OBS_TRANSMISSION_SYSTEMATICS = 0.005
 OBS_OBJECT_TYPE = 'STAR'  # To choose between STAR, HG-AR, MONOCHROMATOR
-OBS_TELESCOPE_TRANSMISSION = 'ctio_throughput.txt'  # telescope transmission file
-OBS_FULL_INSTRUMENT_TRANSMISSON = 'ctio_throughput_300517_v1.txt'  # full instrument transmission file
-OBS_QUANTUM_EFFICIENCY = "qecurve.txt"  # quantum efficiency of the detector file
+OBS_FULL_INSTRUMENT_TRANSMISSON = 'ctio_throughput_300517_v1.txt'  # Full instrument transmission file
+OBS_TRANSMISSION_SYSTEMATICS = 0.005
 OBS_CAMERA_ROTATION = 0  # Camera (x,y) rotation angle with respect to (north-up, east-left) system in degrees
 OBS_CAMERA_DEC_FLIP_SIGN = 1  # Camera (x,y) flip signs with respect to (north-up, east-left) system
 OBS_CAMERA_RA_FLIP_SIGN = 1  # Camera (x,y) flip signs with respect to (north-up, east-left) system

spectractor/simulation/simulator.py

@@ -441,7 +441,7 @@ def simulate(self, A1=1.0, A2=0., A3=0., aerosols=0.05, angstrom_exponent=None,
         -------
         >>> spec = Spectrum("./tests/data/reduc_20170530_134_spectrum.fits")
         >>> spec.disperser.ratio_ratio_order_3over2 = lambda lbda: 0.1
-        >>> psf_poly_params = spec.chromatic_psf.from_table_to_poly_params()
+        >>> psf_poly_params = list(spec.chromatic_psf.from_table_to_poly_params()) * 3
         >>> atmosphere = Atmosphere(airmass=1.2, pressure=800, temperature=10)
         >>> sim = SpectrogramModel(spec, atmosphere=atmosphere, with_background=True, fast_sim=True)
         >>> lambdas, model, model_err = sim.simulate(A2=1, angle=-1.5, psf_poly_params=psf_poly_params)

spectractor/simulation/throughput.py

@@ -93,8 +93,8 @@ def plot_transmission_simple(ax, lambdas, transmissions,  uncertainties=None, la
         >>> plot_transmission_simple(ax, lambdas, transmissions, errors, title="CTIO", label="FGB37")
         >>> lambdas, transmissions, errors = load_transmission(os.path.join(parameters.THROUGHPUT_DIR, "RG715.txt"))
         >>> plot_transmission_simple(ax, lambdas, transmissions, errors, title="CTIO", label="RG715")
-        >>> lambdas, transmissions, errors = load_transmission(os.path.join(parameters.THROUGHPUT_DIR, parameters.OBS_TELESCOPE_TRANSMISSION))
-        >>> plot_transmission_simple(ax, lambdas, transmissions, errors, title="CTIO", label="Telescope")
+        >>> lambdas, transmissions, errors = load_transmission(os.path.join(parameters.THROUGHPUT_DIR, parameters.OBS_FULL_INSTRUMENT_TRANSMISSON))
+        >>> plot_transmission_simple(ax, lambdas, transmissions, errors, title="CTIO", label="Full instrument")
         >>> if parameters.DISPLAY: plt.show()
 
     """