Skip to content

Commit

Permalink
Merge pull request #141 from LSSTDESC/v3.0.2
Browse files Browse the repository at this point in the history
update the timeouts and repair docs
  • Loading branch information
jeremyneveu authored Nov 1, 2023
2 parents 625b659 + eabfc2f commit 75fd02b
Show file tree
Hide file tree
Showing 7 changed files with 46 additions and 38 deletions.
4 changes: 2 additions & 2 deletions config/auxtel.ini
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
3 changes: 2 additions & 1 deletion docs/requirements.txt
Original file line number Diff line number Diff line change
Expand Up @@ -7,4 +7,5 @@ coloredlogs
recommonmark
mistune==2.0.3
docutils<1.0,>=0.16
sphinx-mdinclude
sphinx-mdinclude
sphinx-rtd-theme
2 changes: 1 addition & 1 deletion setup.py
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
18 changes: 11 additions & 7 deletions spectractor/extractor/chromaticpsf.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
48 changes: 23 additions & 25 deletions spectractor/extractor/psf.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
----------
Expand All @@ -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
----------
Expand Down Expand Up @@ -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
Expand All @@ -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
-------
Expand Down Expand Up @@ -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
----------
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
-------
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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
----------
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
3 changes: 3 additions & 0 deletions spectractor/extractor/spectroscopy.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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()
"""
Expand Down
6 changes: 4 additions & 2 deletions spectractor/extractor/spectrum.py
Original file line number Diff line number Diff line change
Expand Up @@ -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')
Expand All @@ -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')
Expand Down

0 comments on commit 75fd02b

Please sign in to comment.