spectractor.fit package

Submodules

spectractor.fit.fitter module

class spectractor.fit.fitter.FitParameter(value: float, label: str, axis_name: str, bounds: list, fixed: bool, err: float, truth: float | None = None)

Bases: object

Container for a parameter to fit on data with FitWorkspace.

value

Array containing the parameter values.

Type:

float

label

Parameter labels for screen print.

Type:

str

axis_name

Parameter labels for plot print.

Type:

str

bounds

2-element list giving the (lower, upper) bounds for the parameter.

Type:

list

fixed

Boolean indicating whether the parameter is fixed.

Type:

bool

err

Uncertainty on parameter value.

Type:

float

truth

Truth value.

Type:

float, optional

Examples

>>> from spectractor.fit.fitter import FitParameter
>>> p = FitParameter(value=1, label="x", axis_name="$x$", bounds=[0, 1], fixed=False, err=0.1, truth=1)
>>> p
x: 1.0 +/- 0.1 (truth=1)
>>> p = FitParameter(value=0.01234, label="t", axis_name="$t$", bounds=[-1, 1], fixed=False, err=0.02)
>>> p
t: 0.01 +/- 0.02
>>> p = FitParameter(value=1, label="x", axis_name="$x$", bounds=[0, 1], fixed=True, err=0)
>>> p
x: 1 (fixed)

axis_name: str

bounds: list

err: float

fixed: bool

label: str

truth: float | None = None

value: float

class spectractor.fit.fitter.FitParameters(values: ndarray | list, labels: list | None = None, axis_names: list | None = None, bounds: list | None = None, fixed: list | None = None, truth: list | None = None, filename: str | None = '', extra: dict | None = None)

Bases: object

Container for the parameters to fit on data with FitWorkspace.

values

Array containing the parameter values.

Type:

np.ndarray

labels

List of the parameter labels for screen print. If None, make a default list with parameters labelled par (default: None).

Type:

list, optional

axis_names

List of the parameter labels for plot print. If None, make a default list with parameters labelled par (default: None).

Type:

list, optional

bounds

List of 2-element list giving the (lower, upper) bounds for every parameter. If None, make a default list with np.infinity boundaries (default: None).

Type:

list, optional

fixed

List of boolean: True to fix a parameter, False to let it free. If None, make a default list with False values: all parameters are free (default: None).

Type:

list, optional

truth

Array of truth parameters (default: None).

Type:

np.ndarray, optional

filename

File name associated to the fitted parameters (usually _spectrum.fits file name) (default: ‘’).

Type:

str, optional

Examples

>>> from spectractor.fit.fitter import FitParameters
>>> params = FitParameters(values=[1, 1, 1, 1, 1])
>>> params.ndim
5
>>> params.values
array([1., 1., 1., 1., 1.])
>>> params.labels
['par0', 'par1', 'par2', 'par3', 'par4']
>>> params.bounds
[[-inf, inf], [-inf, inf], [-inf, inf], [-inf, inf], [-inf, inf]]

axis_names: list | None = None

bounds: list | None = None

property err

Uncertainties on fitted parameters, as the square root of the covariance matrix diagonal.

Returns:

err – The uncertainty array.

Return type:

np.ndarray

Examples

>>> from spectractor.fit.fitter import FitParameters
>>> import numpy as np
>>> params = FitParameters(values=[1, 1, 1, 1, 1], fixed=[False, True, False, True, False])
>>> params.cov = np.array([[1,-0.5,0],[-0.5,4,-1],[0,-1,9]])
>>> params.err
array([1., 0., 2., 0., 3.])

extra: dict | None = None

filename: str | None = ''

fixed: list | None = None

get_fixed_parameters()

Return indices array of fixed parameters.

Examples

>>> params = FitParameters(values=[1, 1, 1, 1, 1], fixed=None)
>>> params.fixed
[False, False, False, False, False]
>>> params.get_fixed_parameters()
array([], dtype=int64)
>>> params = FitParameters(values=[1, 1, 1, 1, 1], fixed=[True, False, True, False, True])
>>> params.fixed
[True, False, True, False, True]
>>> params.get_fixed_parameters()
array([0, 2, 4])

get_free_parameters()

Return indices array of free parameters.

Examples

>>> params = FitParameters(values=[1, 1, 1, 1, 1], fixed=None)
>>> params.fixed
[False, False, False, False, False]
>>> params.get_free_parameters()
array([0, 1, 2, 3, 4])
>>> params = FitParameters(values=[1, 1, 1, 1, 1], fixed=[True, False, True, False, True])
>>> params.fixed
[True, False, True, False, True]
>>> params.get_free_parameters()
array([1, 3])

get_index(label)

Get parameter index given its label.

Parameters:

label (str) – The parameter label.

Returns:

index – The parameter index.

Return type:

int

Examples

>>> params = FitParameters(values=[1, 2, 3, 4], labels=["x", "y", "z", "t"], fixed=[True, False, True, False])
>>> params.get_index("z")
2

get_parameter(key)

Return a FitParameter instance. key argument can be the parameter label or its index value.

Parameters:

key (str, int) – The parameter key.

Returns:

param – A FitParameter instance containing all information about a parameter.

Return type:

FitParameter

Examples

>>> params = FitParameters(values=[1, 2, 3, 4], labels=["x", "y", "z", "t"], fixed=[True, False, True, False])
>>> params.cov = np.array([[1, -0.5], [-0.5, 4]])
>>> p3 = params.get_parameter("t")
>>> p3
t: 4 +/- 2
>>> p0 = params.get_parameter(0)
>>> p0
x: 1.0 (fixed)

property json_filename

labels: list | None = None

property ndim

Number of parameters.

Returns:

ndim

Return type:

int

Examples

>>> from spectractor.fit.fitter import FitParameters
>>> params = FitParameters(values=[1, 1, 1, 1, 1])
>>> params.ndim
5

property nfixed

Number of fixed parameters.

Returns:

nfixed

Return type:

int

Examples

>>> from spectractor.fit.fitter import FitParameters
>>> params = FitParameters(values=[1, 1, 1, 1, 1], fixed=[True, False, True, False, True])
>>> params.nfixed
3

property nfree

Number of free parameters.

Returns:

nfree

Return type:

int

Examples

>>> from spectractor.fit.fitter import FitParameters
>>> params = FitParameters(values=[1, 1, 1, 1, 1], fixed=[True, False, True, False, True])
>>> params.nfree
2

plot_correlation_matrix(live_fit=False)

Compute and plot a correlation matrix.

Save the plot if parameters.SAVE is True. The output file name is build from self.file_name, adding the suffix _correlation.pdf.

Parameters:

live_fit (bool, optional, optional) – If True, model, data and residuals plots are made along the fitting procedure (default: False).

Examples

>>> from spectractor.fit.fitter import FitParameters
>>> params = FitParameters(values=[1, 1, 1], axis_names=["x", "y", "z"])
>>> params.cov = np.array([[1,-0.5,0],[-0.5,1,-1],[0,-1,1]])
>>> params.plot_correlation_matrix()

print_parameters_summary()

Build a string with the best fitting parameters to display on screen. Labels are from self.labels.

Returns:

txt – The printed text.

Return type:

str

Examples

>>> parameters.VERBOSE = True
>>> params = FitParameters(values=[1, 2, 3, 4], labels=["x", "y", "z", "t"], fixed=[True, False, True, False])
>>> params.cov = np.array([[1, -0.5], [-0.5, 4]])
>>> _ = params.print_parameters_summary()

property rho

Correlation matrix computed from the covariance matrix

Returns:

rho – The correlation matrix array.

Return type:

np.ndarray

Examples

>>> from spectractor.fit.fitter import FitParameters
>>> params = FitParameters(values=[1, 1, 1], axis_names=["x", "y", "z"])
>>> params.cov = np.array([[2,-0.5,0],[-0.5,2,-1],[0,-1,2]])
>>> params.rho
array([[ 1.  , -0.25,  0.  ],
       [-0.25,  1.  , -0.5 ],
       [ 0.  , -0.5 ,  1.  ]])

set(label, value)

Set value parameter given its label.

Parameters:

• label (str) – The parameter label.

• value (float) – The new parameter value.

Examples

>>> params = FitParameters(values=[1, 2, 3, 4], labels=["x", "y", "z", "t"], fixed=[True, False, True, False])
>>> params["z"]
3.0
>>> params.set("z", 0)
>>> params["z"]
0.0

truth: list | None = None

property txt_filename

values: ndarray | list

write_json()

write_text(header='')

Save the best fitting parameter summary in a text file.

The file name is build from self.file_name, adding the suffix _bestfit.txt.

Parameters:

header (str, optional) – A header to add to the file (default: “”).

Examples

>>> params = FitParameters(values=[1, 2, 3, 4], labels=["x", "y", "z", "t"],  fixed=[True, False, True, False], filename="test_spectrum.fits")
>>> params.cov = np.array([[1,-0.5,0],[-0.5,1,-1],[0,-1,1]])
>>> params.write_text(header="chi2: 1")

class spectractor.fit.fitter.FitWorkspace(params=None, data=None, err=None, data_cov=None, epsilon=None, file_name='', verbose=False, plot=False, live_fit=False, truth=None)

Bases: object

chisq(p, model_output=False)

Compute the chi square for a set of model parameters p.

Four cases are implemented: diagonal W, 2D W, array of diagonal Ws, array of 2D Ws. The two latter cases are for multiple independent data vectors with W being block diagonal.

Parameters:

• p (array_like) – The array of model parameters.

• model_output (bool, optional) – If true, the simulated model is output.

Returns:

chisq – The chi square value.

Return type:

float

compute_W_with_model_error(model_err)

Propagate model uncertainties to weight matrix W. The method add the mode uncertainties in quadrature to the inverse of the weight matrix W self.data_cov and re-invert it.

Parameters:

model_err (np.ndarray) – Flat array of model uncertainties.

Returns:

W – Weight matrix with model uncertainties propagated.

Return type:

array_like

Examples

1D case:

>>> w = FitWorkspace()
>>> w.data_cov = 2 * np.ones(4)
>>> w.prepare_weight_matrices()
>>> w.compute_W_with_model_error(np.sqrt(2) * np.ones(4))
array([0.25, 0.25, 0.25, 0.25])

2D case:

>>> w = FitWorkspace()
>>> w.data = np.array([1,2,3])
>>> w.data_cov = np.diag([1,1,1])
>>> w.prepare_weight_matrices()
>>> w.compute_W_with_model_error(np.sqrt(3) * np.ones(3))
array([[0.25, 0.  , 0.  ],
       [0.  , 0.25, 0.  ],
       [0.  , 0.  , 0.25]])

Add outliers:

>>> w.outliers = np.array([2])
>>> w.prepare_weight_matrices()
>>> w.compute_W_with_model_error(np.sqrt(3) * np.ones(3))
array([[0.25, 0.  , 0.  ],
       [0.  , 0.25, 0.  ],
       [0.  , 0.  , 0.  ]])

Use sparse matrix:

>>> w.W = scipy.sparse.diags(np.ones(3)).tocsr()
>>> w.W[-1, -1] = 0  # mask last data point
>>> w.W.toarray()
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 0.]])
>>> w.compute_W_with_model_error(np.sqrt(3) * np.ones(3)).toarray()
array([[0.25, 0.  , 0.  ],
       [0.  , 0.25, 0.  ],
       [0.  , 0.  , 0.  ]])

get_bad_indices()

List of indices that are outliers rejected by a sigma-clipping method or other masking method.

Returns:

outliers

Return type:

list

Examples

>>> w = FitWorkspace()
>>> w.data = np.array([np.array([1,2,3]), np.array([1,2,3,4])], dtype=object)
>>> w.outliers = [2, 6]
>>> w.get_bad_indices()
[array([2]), array([3])]

hessian(params)

Experimental function to compute the hessian of a model.

Parameters:

params (array_like) – The array of model parameters.

jacobian(params, model_input=None)

Generic function to compute the Jacobian matrix of a model, with numerical derivatives.

Parameters:

• params (array_like) – The array of model parameters.

• model_input (array_like, optional) – A model input as a list with (x, model, model_err) to avoid an additional call to simulate().

Returns:

J – The Jacobian matrix.

Return type:

np.array

lnlike(p)

Compute the logarithmic likelihood for a set of model parameters p as -0.5*chisq.

Parameters:

p (array_like) – The array of model parameters.

Returns:

lnlike – The logarithmic likelihood value.

Return type:

float

plot_fit()

Generic function to plot the result of the fit for 1D curves.

Returns:

fig – The figure.

Return type:

plt.FigureClass

plot_gradient_descent()

plot_triangle_fisher(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’.

prepare_weight_matrices()

Compute weight matrix \(\mathbf{W}\) self.W as the inverse of data covariance matrix self.data_cov. Cancel weights of data outliers given by self.outliers.

Examples

1D case:

>>> w = FitWorkspace()
>>> w.data_cov = 2 * np.ones(4)
>>> w.prepare_weight_matrices()
>>> w.W
array([0.5, 0.5, 0.5, 0.5])

2D case:

>>> w = FitWorkspace()
>>> w.data = np.array([1,2,3])
>>> w.data_cov = np.diag([1,2,4])
>>> w.prepare_weight_matrices()
>>> w.W
array([[1.  , 0.  , 0.  ],
       [0.  , 0.5 , 0.  ],
       [0.  , 0.  , 0.25]])

Add outliers:

>>> w.outliers = np.array([2])
>>> w.prepare_weight_matrices()
>>> w.W
array([[1. , 0. , 0. ],
       [0. , 0.5, 0. ],
       [0. , 0. , 0. ]])

Use sparse matrix:

>>> w.W = scipy.sparse.diags(np.ones(3))
>>> w.W.toarray()
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])
>>> w.prepare_weight_matrices()
>>> w.W.getformat()
'dia'
>>> w.W.toarray()
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 0.]])

save_gradient_descent()

set_epsilon()

simulate(*p)

Compute the model prediction given a set of parameters.

Parameters:

p (array_like) – Array of parameters for the computation of the model.

Returns:

• x (array_like) – The abscisse of the model prediction.

• model (array_like) – The model prediction.

• model_err (array_like) – The uncertainty on the model prediction.

Examples

>>> w = FitWorkspace()
>>> p = np.zeros(3)
>>> x, model, model_err = w.simulate(*p)

weighted_residuals(p)

Compute the weighted residuals array for a set of model parameters p.

Parameters:

p (array_like) – The array of model parameters.

Returns:

residuals – The array of weighted residuals.

Return type:

np.array

class spectractor.fit.fitter.RegFitWorkspace(w, opt_reg=0.01, verbose=False, live_fit=False)

Bases: FitWorkspace

plot_fit()

Generic function to plot the result of the fit for 1D curves.

Returns:

fig – The figure.

Return type:

plt.FigureClass

print_regularisation_summary()

run_regularisation(Ndof=None)

set_regularisation_with_Ndof(Ndof)

Find regularisation parameter $r$ that checks $Tr(R) = Ndof$.

Parameters:

Ndof (float) – Number of degrees of freedom, ie $Tr(R)$.

Returns:

r – Regularisation parameter.

Return type:

float

simulate(log10_r)

Compute the model prediction given a set of parameters.

Parameters:

p (array_like) – Array of parameters for the computation of the model.

Returns:

• x (array_like) – The abscisse of the model prediction.

• model (array_like) – The model prediction.

• model_err (array_like) – The uncertainty on the model prediction.

Examples

>>> w = FitWorkspace()
>>> p = np.zeros(3)
>>> x, model, model_err = w.simulate(*p)

spectractor.fit.fitter.gradient_descent(fit_workspace, niter=10, xtol=0.001, ftol=0.001, with_line_search=True)

Four cases are implemented: diagonal W, 2D W, array of diagonal Ws, array of 2D Ws. The two latter cases are for multiple independent data vectors with W being block diagonal.

Parameters:

• fit_workspace (FitWorkspace)

• niter

• xtol

• ftol

• with_line_search

spectractor.fit.fitter.read_fitparameter_json(json_filename)

Read JSON file and store data in FitParameters instance.

Parameters:

json_filename (str) – The JSON file name.

Returns:

params – A FitParameters instance to loaded from JSON json_filename.

Return type:

FitParameters

Examples

>>> params = FitParameters(values=[1, 2, 3, 4], labels=["x", "y", "z", "t"],  fixed=[True, False, True, False], filename="test_spectrum.fits")
>>> params.cov = np.array([[1,-0.5,0],[-0.5,1,-1],[0,-1,1]])
>>> _ = write_fitparameter_json(params.json_filename, params, extra={"chi2": 1})
>>> new_params = read_fitparameter_json(params.json_filename)
>>> new_params.values
array([1., 2., 3., 4.])

spectractor.fit.fitter.run_gradient_descent(fit_workspace, xtol, ftol, niter, verbose=False, with_line_search=True)

spectractor.fit.fitter.run_minimisation(fit_workspace, method='newton', xtol=0.0001, ftol=0.0001, niter=50, verbose=False, with_line_search=True, minimizer_method='L-BFGS-B')

spectractor.fit.fitter.run_minimisation_sigma_clipping(fit_workspace, method='newton', xtol=0.0001, ftol=0.0001, niter=50, sigma_clip=5.0, niter_clip=3, verbose=False, with_line_search=True)

spectractor.fit.fitter.run_simple_newton_minimisation(fit_workspace, xtol=1e-08, ftol=1e-08, niter=50, verbose=False)

spectractor.fit.fitter.simple_newton_minimisation(fit_workspace, niter=10, xtol=0.001, ftol=0.001)

Experimental function to minimize a function.

Parameters:

• fit_workspace (FitWorkspace)

• niter

• xtol

• ftol

spectractor.fit.fitter.write_fitparameter_json(json_filename, params, extra=None)

Save FitParameters attributes as a JSON file.

Parameters:

• json_filename (str) – JSON file name.

• params (FitParameters) – A FitParameters instance to save in JSON json_filename.

• extra (dict, optional) – Extra information to write in the JSON file.

Returns:

jsontxt – The JSON dictionnary as string.

Return type:

str

Examples

>>> params = FitParameters(values=[1, 2, 3, 4], labels=["x", "y", "z", "t"],  fixed=[True, False, True, False], filename="test_spectrum.fits")
>>> params.cov = np.array([[1,-0.5,0],[-0.5,1,-1],[0,-1,1]])
>>> jsonstr = write_fitparameter_json(params.json_filename, params, extra={"chi2": 1})
>>> jsonstr  
'{"values": [1.0, 2.0, 3.0, 4.0], "labels": ["x", "y", "z", "t"],..."extra": {"chi2": 1}...

spectractor.fit.statistics module

spectractor.fit.fit_spectrogram module

class spectractor.fit.fit_spectrogram.SpectrogramFitWorkspace(spectrum, atmgrid_file_name='', fit_angstrom_exponent=False, verbose=False, plot=False, live_fit=False, truth=None)

Bases: FitWorkspace

crop_spectrogram()

Crop the spectrogram in the middle, keeping a vertical width of 2*parameters.PIXWIDTH_SIGNAL around the signal region.

get_spectrogram_truth()

Load the truth parameters (if provided) from the file header.

jacobian(params, model_input=None)

Generic function to compute the Jacobian matrix of a model, with numerical derivatives.

Parameters:

• params (array_like) – The array of model parameters.

• model_input (array_like, optional) – A model input as a list with (x, model, model_err) to avoid an additional call to simulate().

Returns:

J – The Jacobian matrix.

Return type:

np.array

plot_fit()

Plot the fit result.

Examples

>>> spec = Spectrum('tests/data/reduc_20170530_134_spectrum.fits')
>>> parameters.SPECTRACTOR_ATMOSPHERE_SIM = "libradtran"
>>> w = SpectrogramFitWorkspace(spec, verbose=True)
>>> lambdas, model, model_err = w.simulate(*w.params.values)
>>> w.plot_fit()

from spectractor.fit.fit_spectrogram import SpectrogramFitWorkspace
file_name = 'tests/data/reduc_20170530_134_spectrum.fits'
atmgrid_file_name = file_name.replace('spectrum', 'atmsim')
fit_workspace = SpectrogramFitWorkspace(file_name, atmgrid_file_name=atmgrid_file_name, verbose=True)
A1, A2, aerosols, ozone, pwv, D, shift_x, shift_y, angle, *psf = fit_workspace.p
lambdas, model, model_err = fit_workspace.simulation.simulate(A1, A2, aerosols, ozone, pwv, D, shift_x,
                                                              shift_y, angle, psf)
fit_workspace.lambdas = lambdas
fit_workspace.model = model
fit_workspace.model_err = model_err
fit_workspace.plot_fit()

(Source code)

plot_spectrogram_comparison_simple(ax, title='', extent=None, dispersion=False)

Method to plot a spectrogram issued from data and compare it with simulations.

Parameters:

• ax (Axes) – Axes instance of shape (3, 2).

• title (str, optional) – Title for the simulation plot (default: ‘’).

• extent (array_like, optional) – Extent argument for imshow to crop plots (default: None).

• dispersion (bool, optional) – If True, plot a colored bar to see the associated wavelength color along the x axis (default: False).

set_mask(params=None)

Parameters:

params

Examples

>>> spec = Spectrum('tests/data/reduc_20170530_134_spectrum.fits')
>>> parameters.SPECTRACTOR_ATMOSPHERE_SIM = "libradtran"
>>> w = SpectrogramFitWorkspace(spec, verbose=True)
>>> _ = w.simulate(*w.params.values)
>>> w.plot_fit()

simulate(*params)

Interface method to simulate a spectrogram.

Parameters:

params (array_like) – Simulation parameter array.

Returns:

• lambdas (array_like) – Array of wavelengths (1D).

• model (array_like) – Flat 1D array of the spectrogram simulation.

• model_err (array_like) – Flat 1D array of the spectrogram simulation uncertainty.

Examples

>>> spec = Spectrum('tests/data/reduc_20170530_134_spectrum.fits')
>>> parameters.SPECTRACTOR_ATMOSPHERE_SIM = "libradtran"
>>> w = SpectrogramFitWorkspace(spec, verbose=True)
>>> lambdas, model, model_err = w.simulate(*w.params.values)
>>> w.plot_fit()

spectractor.fit.fit_spectrogram.lnprob_spectrogram(p)

Logarithmic likelihood function to maximize in MCMC exploration.

Parameters:

p (array_like) – Array of SpectrogramFitWorkspace parameters.

Returns:

lp – Log of the likelihood function.

Return type:

float

spectractor.fit.fit_spectrogram.run_spectrogram_minimisation(fit_workspace, method='newton', verbose=False)

Interface function to fit spectrogram simulation parameters to data.

Parameters:

• fit_workspace (SpectrogramFitWorkspace) – An instance of the SpectrogramFitWorkspace class.

• method (str, optional) – Fitting method (default: ‘newton’).

Examples

>>> spec = Spectrum('tests/data/reduc_20170530_134_spectrum.fits')
>>> parameters.SPECTRACTOR_ATMOSPHERE_SIM = "libradtran"
>>> w = SpectrogramFitWorkspace(spec, verbose=True, atmgrid_file_name='tests/data/reduc_20170530_134_atmsim.fits')
>>> parameters.VERBOSE = True
>>> run_spectrogram_minimisation(w, method="newton")

Module contents