diff --git a/config/auxtel.ini b/config/auxtel.ini index 3e791d57b..babf21ccc 100644 --- a/config/auxtel.ini +++ b/config/auxtel.ini @@ -16,9 +16,9 @@ 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 +SPECTRACTOR_FIT_TIMEOUT_PER_ITER = 1200 # maximum time per gradient descent before TimeoutError in seconds -SPECTRACTOR_FIT_TIMEOUT = 3600 +SPECTRACTOR_FIT_TIMEOUT = 7200 [instrument] # instrument name diff --git a/docs/requirements.txt b/docs/requirements.txt index 339504629..553dfd0f7 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -7,4 +7,5 @@ coloredlogs recommonmark mistune==2.0.3 docutils<1.0,>=0.16 -sphinx-mdinclude \ No newline at end of file +sphinx-mdinclude +sphinx-rtd-theme \ No newline at end of file diff --git a/setup.py b/setup.py index 0e992b5ff..cc9386411 100644 --- a/setup.py +++ b/setup.py @@ -8,7 +8,7 @@ if os.getenv('CONDA_PREFIX') is None: reqs = open('requirements.txt', 'r').read().strip().splitlines() -if os.getenv('READTHEDOCS'): +if os.getenv('READTHEDOCS') and 'mpi4py' in reqs: reqs.remove('mpi4py') with open('README.md') as file: diff --git a/spectractor/extractor/chromaticpsf.py b/spectractor/extractor/chromaticpsf.py index 6c3e3ce80..832ec833c 100644 --- a/spectractor/extractor/chromaticpsf.py +++ b/spectractor/extractor/chromaticpsf.py @@ -2406,32 +2406,36 @@ def amplitude_derivatives(self): def jacobian(self, params, epsilon, model_input=None): r"""Generic function to compute the Jacobian matrix of a model, linear parameters being fixed (see Notes), - with analytical or numerical derivatives. Analytical derivatives are performed if self.analytical is True. - Let's write :math:`\theta` the non-linear model parameters. If the model is written as: + with analytical or numerical derivatives. Analytical derivatives are performed if `self.analytical` is True. + Let's write :math:`\theta` the non-linear model parameters. If the model is written as: + .. math:: \mathbf{I} = \mathbf{M}(\theta) \cdot \hat{\mathbf{A}}(\theta), this jacobian function returns: + .. math:: \frac{\partial \mathbf{I}}{\partial \theta} = \frac{\partial \mathbf{M}}{\partial \theta} \cdot \hat{\mathbf{A}}. Notes ----- - The gradient descent is performed on the non-linear parameters :math:`\theta` (PSF shape and position). Linear parameters :math:`\mathbf{A}`(amplitudes) are computed on the fly. + The gradient descent is performed on the non-linear parameters :math:`\theta` (PSF shape and position). Linear parameters :math:`\mathbf{A}` (amplitudes) are computed on the fly. Therefore, :math:`\chi^2` is a function of :math:`\theta` only + .. math :: \chi^2(\theta) = \chi'^2(\theta, \hat{\mathbf{A}} - whose partial derivatives on :math:`\theta` for gradient descent are + whose partial derivatives on :math:`\theta` for gradient descent are: + .. math :: - \frac{\partial \chi^2}{\partial \theta} = \left.\left(\frac{\partial \chi'^2}{\partial \theta} + \frac{\partial \chi'^2}{\partial \mathbf{A}} \frac{\partial \mathbf{A}}{\partial \theta}\right)\right\vert_{\mathbf{A} = \hat \mathbf{A}} + \frac{\partial \chi^2}{\partial \theta} = \left.\left(\frac{\partial \chi'^2}{\partial \theta} + \frac{\partial \chi'^2}{\partial \mathbf{A}} \frac{\partial \mathbf{A}}{\partial \theta}\right)\right\vert_{\mathbf{A} = \hat{\mathbf{A}}} - By definition, :math:`\left.\partial \chi'^2/\partial \mathbf{A}\right\vert_{\mathbf{A} = \hat \mathbf{A}}=0` then :math:`\chi^2` partial derivatives must be performed with fixed :math:`\mathbf{A} = \hat \mathbf{A}` - for gradient descent. `self.amplitude_priors_method` is temporarily switched to "keep" in `self.jacobian()` to use previously computed :math:`\hat \mathbf{A}` solution. + By definition, :math:`\left.\partial \chi'^2/\partial \mathbf{A}\right\vert_{\mathbf{A} = \hat{\mathbf{A}}}=0` then :math:`\chi^2` partial derivatives must be performed with fixed :math:`\mathbf{A} = \hat{\mathbf{A}}` + for gradient descent. `self.amplitude_priors_method` is temporarily switched to "keep" in `self.jacobian()` to use previously computed :math:`\hat{\mathbf{A}}` solution. Parameters diff --git a/spectractor/extractor/psf.py b/spectractor/extractor/psf.py index 020f03ee7..63acf2bc4 100644 --- a/spectractor/extractor/psf.py +++ b/spectractor/extractor/psf.py @@ -18,10 +18,10 @@ def evaluate_moffat1d_normalisation(gamma, alpha): .. math :: - A = \frac{\Gamma(alpha)}{\gamma \sqrt{\pi} \Gamma(alpha -1/2)} + A = \frac{\Gamma(\alpha)}{\gamma \sqrt{\pi} \Gamma(\alpha -1/2)} \quad\text{with}, \alpha > 1/2 - Note that this function is defined only for :math:`alpha > 1/2`. + Note that this function is defined only for :math:`\alpha > 1/2`. Parameters ---------- @@ -48,10 +48,10 @@ def evaluate_moffat1d_normalisation_dalpha(norm, alpha): .. math :: - A = \frac{\Gamma(alpha)}{\gamma \sqrt{\pi} \Gamma(alpha -1/2)} + A = \frac{\Gamma(\alpha)}{\gamma \sqrt{\pi} \Gamma(\alpha -1/2)} \quad\text{with}, \alpha > 1/2 - Note that this function is defined only for :math:`alpha > 1/2`. + Note that this function is defined only for :math:`\alpha > 1/2`. Parameters ---------- @@ -83,8 +83,8 @@ def evaluate_moffat1d(y, amplitude, y_c, gamma, alpha, norm): # pragma: no cove f(y) \propto \frac{A}{\left[ 1 +\left(\frac{y-y_c}{\gamma}\right)^2 \right]^\alpha} \quad\text{with}, \alpha > 1/2 - Note that this function is defined only for :math:`alpha > 1/2`. The normalisation factor - :math:`\frac{\Gamma(alpha)}{\gamma \sqrt{\pi} \Gamma(alpha -1/2)}` is not included as special functions + Note that this function is defined only for :math:`\alpha > 1/2`. The normalisation factor + :math:`\frac{\Gamma(\alpha)}{\gamma \sqrt{\pi} \Gamma(\alpha -1/2)}` is not included as special functions are not supported by numba library. Parameters @@ -100,7 +100,7 @@ def evaluate_moffat1d(y, amplitude, y_c, gamma, alpha, norm): # pragma: no cove alpha: float Exponent :math:`\alpha` of the Moffat function. norm: float - Normalisation :math:`\frac{\Gamma(alpha)}{\gamma \sqrt{\pi} \Gamma(alpha -1/2)}`. + Normalisation :math:`\frac{\Gamma(\alpha)}{\gamma \sqrt{\pi} \Gamma(\alpha -1/2)}`. Returns ------- @@ -156,10 +156,10 @@ def evaluate_moffat1d_jacobian(y, amplitude, y_c, gamma, alpha, norm, dnormda, f .. math :: - f(y) \propto \frac{A}{\left[ 1 +\left(\frac{y-y_c}{\gamma}\right)^2 \right]^\alpha}\times \frac{\Gamma(alpha)}{\gamma \sqrt{\pi} \Gamma(alpha -1/2)} + f(y) \propto \frac{A}{\left[ 1 +\left(\frac{y-y_c}{\gamma}\right)^2 \right]^\alpha}\times \frac{\Gamma(\alpha)}{\gamma \sqrt{\pi} \Gamma(\alpha -1/2)} \quad\text{with}, \alpha > 1/2 - Note that this function is defined only for :math:`alpha > 1/2`. + Note that this function is defined only for :math:`\alpha > 1/2`. Parameters ---------- @@ -174,7 +174,7 @@ def evaluate_moffat1d_jacobian(y, amplitude, y_c, gamma, alpha, norm, dnormda, f alpha: float Exponent :math:`\alpha` of the Moffat function. norm: float - Normalisation :math:`\frac{\Gamma(alpha)}{\gamma \sqrt{\pi} \Gamma(alpha -1/2)}`. + Normalisation :math:`\frac{\Gamma(\alpha)}{\gamma \sqrt{\pi} \Gamma(\alpha -1/2)}`. dnormda: float Derivatives of the normalisation with respect to alpha. fixed: array_like @@ -241,8 +241,8 @@ def evaluate_moffatgauss1d(y, amplitude, y_c, gamma, alpha, eta_gauss, sigma, no \frac{1}{\left[ 1 +\left(\frac{y-y_c}{\gamma}\right)^2 \right]^\alpha}+ \eta e^{-(y-y_c)^2/(2\sigma^2)}\right\rbrace \quad\text{ and } \quad \eta < 0, \alpha > 1/2 - Note that this function is defined only for :math:`alpha > 1/2`. The normalisation factor for the Moffat+Gauss - :math:`\frac{\Gamma(alpha)}{\gamma \sqrt{\pi} \Gamma(alpha -1/2)} + \eta \sqrt{2\pi} \sigma` is not included as special functions + Note that this function is defined only for :math:`\alpha > 1/2`. The normalisation factor for the Moffat+Gauss + :math:`\frac{\Gamma(\alpha)}{\gamma \sqrt{\pi} \Gamma(\alpha -1/2)} + \eta \sqrt{2\pi} \sigma` is not included as special functions are not supproted by the numba library. Parameters @@ -262,7 +262,7 @@ def evaluate_moffatgauss1d(y, amplitude, y_c, gamma, alpha, eta_gauss, sigma, no sigma: float Width :math:`\sigma` of the Gaussian function. norm_moffat: float - Normalisation :math:`\frac{\Gamma(alpha)}{\gamma \sqrt{\pi} \Gamma(alpha -1/2)}`. + Normalisation :math:`\frac{\Gamma(\alpha)}{\gamma \sqrt{\pi} \Gamma(\alpha -1/2)}`. Returns ------- @@ -326,12 +326,12 @@ def evaluate_moffatgauss1d_jacobian(y, amplitude, y_c, gamma, alpha, eta_gauss, .. math :: f(y) \propto A \left\lbrace - \frac{1}{\left[ 1 +\left(\frac{y-y_c}{\gamma}\right)^2 \right]^\alpha} \times \frac{\Gamma(alpha)}{\gamma \sqrt{\pi} \Gamma(alpha -1/2)} + \frac{1}{\left[ 1 +\left(\frac{y-y_c}{\gamma}\right)^2 \right]^\alpha} \times \frac{\Gamma(\alpha)}{\gamma \sqrt{\pi} \Gamma(\alpha -1/2)} - \eta e^{-(y-y_c)^2/(2\sigma^2)}\right\rbrace \quad\text{ and } \quad \eta < 0, \alpha > 1/2 - Note that this function is defined only for :math:`alpha > 1/2`. The normalisation factor for the Moffat - :math:`\frac{\Gamma(alpha)}{\gamma \sqrt{\pi} \Gamma(alpha -1/2)} + \eta \sqrt{2\pi} \sigma` is not included as special functions + Note that this function is defined only for :math:`\alpha > 1/2`. The normalisation factor for the Moffat + :math:`\frac{\Gamma(\alpha)}{\gamma \sqrt{\pi} \Gamma(\alpha -1/2)} + \eta \sqrt{2\pi} \sigma` is not included as special functions are not supproted by the numba library. Parameters @@ -351,7 +351,7 @@ def evaluate_moffatgauss1d_jacobian(y, amplitude, y_c, gamma, alpha, eta_gauss, sigma: float Width :math:`\sigma` of the Gaussian function. norm_moffat: float - Normalisation :math:`\frac{\Gamma(alpha)}{\gamma \sqrt{\pi} \Gamma(alpha -1/2)}`. + Normalisation :math:`\frac{\Gamma(\alpha)}{\gamma \sqrt{\pi} \Gamma(\alpha -1/2)}`. dnormda: float Derivatives of the normalisation with respect to alpha. fixed: array_like @@ -430,7 +430,7 @@ def evaluate_moffat2d(x, y, amplitude, x_c, y_c, gamma, alpha): # pragma: no co \quad\text{with}\quad \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x, y) \mathrm{d}x \mathrm{d}y = A - Note that this function is defined only for :math:`alpha > 1`. + Note that this function is defined only for :math:`\alpha > 1`. Note that the normalisation of a 2D Moffat function is analytical so it is not expected that the sum of the output array is equal to :math:`A`, but lower. @@ -589,7 +589,7 @@ def evaluate_moffatgauss2d(x, y, amplitude, x_c, y_c, gamma, alpha, eta_gauss, s \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x, y) \mathrm{d}x \mathrm{d}y = A \quad\text{and} \quad \eta < 0 - Note that this function is defined only for :math:`alpha > 1`. + Note that this function is defined only for :math:`\alpha > 1`. Parameters ---------- @@ -672,8 +672,7 @@ def evaluate_gauss1d(y, amplitude, y_c, sigma): # pragma: no cover .. math :: - f(x, y) = \frac{A}{\sigma \sqrt{2 \pi}\left\lbrace e^{-\left[ \left(x-x_c\right)^2\right]/(2 \sigma^2)} - \right\rbrace + f(y) = \frac{A}{\sigma \sqrt{2 \pi}\left\lbrace e^{-\left[ \left(y-y_c\right)^2\right]/(2 \sigma^2)}\right\rbrace .. math :: \quad\text{with}\quad @@ -743,8 +742,7 @@ def evaluate_gauss1d_jacobian(y, amplitude, y_c, sigma, fixed): # pragma: no co .. math :: - f(x, y) = \frac{A}{\sigma \sqrt{2 \pi}\left\lbrace e^{-\left[ \left(x-x_c\right)^2\right]/(2 \sigma^2)} - \right\rbrace + f(y) = \frac{A}{\sigma \sqrt{2 \pi}\left\lbrace e^{-\left[ \left(y-y_c\right)^2\right]/(2 \sigma^2)}\right\rbrace .. math :: \quad\text{with}\quad @@ -1206,7 +1204,7 @@ def evaluate(self, pixels, values=None): .. math:: - f(y) \propto \frac{A \Gamma(alpha)}{\gamma \sqrt{\pi} \Gamma(alpha -1/2)}, + f(y) \propto \frac{A \Gamma(\alpha)}{\gamma \sqrt{\pi} \Gamma(\alpha -1/2)}, \quad \int_{y_{\text{min}}}^{y_{\text{max}}} f(y) \mathrm{d}y = A Parameters @@ -1576,7 +1574,7 @@ def evaluate(self, pixels, values=None): .. math:: - f(y) \propto \frac{A}{ \frac{\Gamma(alpha)}{\gamma \sqrt{\pi} \Gamma(alpha -1/2)}+\eta\sqrt{2\pi}\sigma}, + f(y) \propto \frac{A}{ \frac{\Gamma(\alpha)}{\gamma \sqrt{\pi} \Gamma(\alpha -1/2)}+\eta\sqrt{2\pi}\sigma}, \quad \int_{y_{\text{min}}}^{y_{\text{max}}} f(y) \mathrm{d}y = A Parameters diff --git a/spectractor/extractor/spectroscopy.py b/spectractor/extractor/spectroscopy.py index dc774204e..c741b0eb1 100644 --- a/spectractor/extractor/spectroscopy.py +++ b/spectractor/extractor/spectroscopy.py @@ -430,6 +430,7 @@ def build_detected_line_table(self, amplitude_units="", calibration_only=False): -------- Creation of a mock spectrum with emission and absorption lines + >>> from spectractor.extractor.spectrum import Spectrum, detect_lines >>> lambdas = np.arange(300,1000,1) >>> spectrum = 1e4*np.exp(-((lambdas-600)/200)**2) @@ -444,12 +445,14 @@ def build_detected_line_table(self, amplitude_units="", calibration_only=False): >>> fwhm_func = interp1d(lambdas, 0.01 * lambdas) Detect the lines + >>> lines = Lines([HALPHA, HBETA, O2_1], hydrogen_only=True, ... atmospheric_lines=True, redshift=0, emission_spectrum=True) >>> global_chisq = detect_lines(lines, lambdas, spectrum, spectrum_err, fwhm_func=fwhm_func) >>> assert(global_chisq < 1) Print the result + >>> spec.lines = lines >>> t = lines.build_detected_line_table() """ diff --git a/spectractor/extractor/spectrum.py b/spectractor/extractor/spectrum.py index a75cfd12c..46b9fe3a5 100644 --- a/spectractor/extractor/spectrum.py +++ b/spectractor/extractor/spectrum.py @@ -814,7 +814,8 @@ def load_spectrum(self, input_file_name, spectrogram_file_name_override=None, Examples -------- - # Latest Spectractor output format: do not provide a config file (parameters are loaded from file header) + Latest Spectractor output format: do not provide a config file (parameters are loaded from file header) + >>> from spectractor import parameters >>> s = Spectrum(config="") >>> s.load_spectrum('tests/data/reduc_20170530_134_spectrum.fits') @@ -826,7 +827,8 @@ def load_spectrum(self, input_file_name, spectrogram_file_name_override=None, >>> assert parameters.CCD_REBIN == s.header["REBIN"] >>> assert s.parallactic_angle == s.header["PARANGLE"] - # Spectractor output format older than version <=2.3: must give the config file + Spectractor output format older than version <=2.3: must give the config file + >>> parameters.VERBOSE = False >>> s = Spectrum(config="./config/ctio.ini") >>> s.load_spectrum('tests/data/reduc_20170605_028_spectrum.fits')