nDspec API

Docstrings for every class and function in the nDspec modelling software

Operator Class

class ndspec.Operator.nDspecOperator

Generic class from which all nDspec operators are inherited. It contains standardized, private methods to manipulate arrays by performing e.g. interpolation, integration over a range or axis, etc. These can be called by each specific operator to perform physical calculations, such as integrating over a Fourier frequency range to calculate a lag spectrum. This class is not to be instantiated by itself.

Response Matrix Class

class ndspec.Response.ResponseMatrix(resp_path, arf_path=None)

This class handles folding an energy-dependent, multi-dimensional product, through a given X-ray instrument response matrix, including the effects of both the redistribution matrix and the effective area.

Currently the X-ray observatories explicitely supported are NICER, NuSTAR and RXTE. Swift/XRT/XMM/Chandra/Insight-HXMT/Astrosat should be compatible but have not yet been tested. Support for these missions is pending.

XRISM and Athena are NOT compatible with this code due to a) the large size of their responses and b) the unusual format of the contents of the responses when read with Astropy. Fixes for both of these issues are pending.

Parameters:

resp_path: string

The path to the response file (either .rmf or .rsp) to load.

Other parameters:

arf_path: string

Optional path to an effective area file (.arf) to load along with the redistribution matrix. Note that some mission tools produce .rmf files which also contain the telescope effective area, in which case loading the arf is not necessary.

Attributes:

chans: numpy.array(int)

An array of integers of size (n_chans) representing each channel in the response.

n_chans: int

The length of the chans array.

emin, emax: numpy.array(float)

The arrays with the minimum/maximum energies in each channel in the chans array. After rebinning on a new grid, it contains the minimum and maximum energies of each bin in the new grid.

energ_lo, energ_hi: numpy.array(float)

The arrays with the minimum/maximum energies bins that sample the instrument rmf/arf.

n_energs: int

The length of the energ_lo/energ_hi arrays.

resp_matrix: numpy.array(float x float)

The instrument response in matrix form, with size (n_energs x n_chans). Only contains the redistribution matrix if no arf is loaded yet.

specresp: numpy.array(float)

The instrument effective area as a function of energy as contained in the arf file.

energ_rebin: bool

A flag to check whether the response was re-binned over energies, in which case folding requires a renormalization constant to conserve the total photon counts

has_arf: bool

A flag to check whether only a rmf file has been loaded, or a full rmf+arf response. This varies between observatories.

mission: str

A string which tracks the observatory for which the response is defined.

instrument: str

A string which tracks the instrument for which the response is defined.

convolve_response(model_input, units_in='xspec', units_out='kev')

This method applies the response matrix loaded in the class to a user defined mode. Two input model normalizations are supported: “rate” normalization, which assumes the input is in units of count rate, or “xspec” normalization, which assumes the model is in units of count rate times energy bin width. Two output normalizations are supported: “kev”, in which case the model is returned in units of counts/s/kev, or “channel”, in which case the model is returned in units of counts/s/channel.

Parameters:

model_input: np.array(float,float) or CrossSpectrum

Either a) a 2-d array of size (n_energs x arbirtrary length), containing the input model as a function of energy and optionally an additional quantity (Fourier frequency, time, pulse phase, etc.), or b) a CrossSpectrum object from nDspec, containing the model cross spectrum to be folded with the instrument response.

units_in: string, default=”xspec”

A string detailing the normalization of the iinput model. The base “xspec” normalization assumes the input is in units of integrated counts per unit time in each energy bin; “rate” normalization assumes the input is in units of integrated count rate per unit time and energy in each bin.

units_out: string, default=”kev”

A string setting the normalization of the output model. The default “kev” normalizations returns a model in units of counts/s/kev, “channel” instead returns a model in units of counts/s/channel.

Returns:

conv_model, np.array(float,float) or CrossSpectrum

Either a) a 2-d array of size (n_chans x arbitrary length), containing the input model as a function of energy channel and a secondary quantity identical to the input model_input (Fourier frequency, time, pulse phase, etc.), or a CrossSpectrum object from nDspec, containing the folded model cross spectrum as a function of energy channel and Fourier frequency.

diagonal_matrix(num)

Returns a diagonal identity matrix, which by definition contains only ones on the diagonal and zeroes otherwise.

Parameters:

num: int

The dimension of the desired matrix.

Returns:

diag_resp: np.array(float,float)

An identity matrix of size (num x num).

load_arf(filepath)

This method reads an effective area .arf file, and applies it to a redistribution matrix previously loaded with the load_rmf method. The arf-corrected matrix over-writes the class attribute resp_matrix.

Additionally, the array of effective area vs energy is stored in the specresp class attribute.

Parameters:

filepath: string

The path to the .arf file to be loaded.

load_rmf(filepath)

This method loads either a rmf or rsp file (containing the redistribution matrix or full response, respectively), and sets the class attributes from it using astropy.

Parameters:

filepath: string

The path to the .rmf or .rsp file to be loaded.

plot_arf(plot_scale='log', return_plot=False)

Plots the instrument effective area, if one has been loaded, as a function of energy.

Parameters:

plot_scale: string, default=”log”

Switches between log10(arf) (plot_scale=”log”, the default behavior) and just the arf (plot_scale=”lin”).

Returns:

fig: matplotlib.figure, optional

The plot object produced by the method.

plot_response(plot_type='channel', return_plot=False)

Plots the instrument response as a function of incoming energy and instrument channel. For ease of visualization, the z-axis plots the base-10 logarithm of the response matrix.

Parameters:

plot_type: string, default=”channel”

Sets the units of the X-axis to be either the channel number (by default) or the bounds of each channel (plot_type=”energy”).

Returns:

fig: matplotlib.figure, optional

The plot object produced by the method.

rebin_channels(new_bounds_lo, new_bounds_hi)

This method rebins the response matrix resp_matrix to an input, continuous grid of channels, and updates the emin, emax and n_chans attributes appropriately. In order to avoid potential issues with gain shifting due bin edges, the method also re-aligns the input grid such that its bounds match the (finer) existing grid. With this method, both the probability to assign an energy channel to a recorded photon, and the total number of photons in the model, are conserved after rebinning.

Parameters:

new_bounds_lo: np.array(float)

An array of energies with the lower bound of each energy channel.

new_bounds_hi: np.array(float)

An array of energies with the upper bound of each energy channel.

Returns:

bin_resp: ResponseMatrix

A ResponseMatrix object containing the same response loaded in the self object, but rebinned over the channel axis to the input grid.

rebin_energies(factor)

This method rebins the response matrix resp_matrix in energy by grouping the a numer of existing binning (defined by the energ_lo and energ_hi attributes) in coarser bins. The input parameter “factor” controls how many bins of the old grid will be grouped in the new grid. With this method, both the probability to assign an energy channel to a recorded photon, and the total number of photons in the model, are conserved after rebinning. Note that rebinning in energy is VERY dangerous and should only be done very carefully.

Parameters:

factor: integer

The number of bins of the old energy grid that will be grouped in the new energy grid

Returns:

bin_resp: ResponseMatrix

A ResponseMatrix object containing the same response loaded in the self object, but rebinned over the energy axis to the input grid.

set_exposure_time(time)

Adjusts the response matrix to a new time. Some mission’s FITS files do not contain this information and thus need to be set manually.

Parameters

timefloat

Exposure time of observation.

Returns

None.

unfold_response(array, units_in='kev')

Unfolds an array through the instrument response. Note that plotting data in this fashion can be EXTREMELY misleading and should be done with care. In nDspec we define an unfolded model as: unfolded(H) = counts(H)/exposure*sum(rmf*arf), where H is a given bounds in energy channels, exposure is the exposure time of the observation, rmf*arf is the instrument response, and the sum is carried out every the energy bins of the response. Users also need to specify whether the input array is in units of counts/s/keV or counts/s/channel; if they do so correctly, the output of this method is in photon density - counts/s/keV/cm^2. Assuming 0 background, this definition of unfolding is identical to Isis, regardless of model choice, and Xspec, as long as the model is a constant in each energy bin. Converting to flux units - ie, energy/s/area, then requires multiplying the output of this method by H^2, identically to the “eeunfold” method in Xspec. Unlike Isis and Xspec, this method also supports unfolding two-d arrays, e.g. cross spectra.

Parameters:

array: np.array(int,int)

The input array to be unfolded, of size (n_energs,arbitrary). In the x-axis it needs to be defined over the instrument energy grid, in the y-axis it can be any size (e.g., over a grid of Fourier frequencies).

Returns:

unfold_model: np.array(float,float)

The array unfolded through the instrument response, of size (n_chans,arbitrary), defined over the energy bounds of each channel in the response. The y-axis is identical to the input.

Timing Classes

class ndspec.Timing.FourierProduct(times, freqs=0, method='fft', verbose=False)

Parent class for all Fourier operators in the software. It is mainly used to handle performing Fourier transforms from the time to the frequency domain, which is implemented using two methods. It is not to be instantiated by itself.

The first simply uses the numpy fftw to calculate a standard fast Fourier transform; the second uses the sinc function decomposition described in Uttley and Malzac 2023 (arXiv:2312.08302), eq. 27.

Depending on the input provided to the constructor, the class can switch automatically between the two implementations to ensure the Fourier transform calculation is correct.

Parameters:

times: np.array(float)

The array of times over which the quantity to be transformed is defined. If using the fft method, this array also defines the frequencies over which the transform is computed.

freqs: np.array(float)

Defines the range of Fourier frequencies over which the input quantity is to be transformed, only if using the sinc function method. Has no effect with fft.

method: string{ “fft” | “sinc” }, default “fft”

The computational method to calculate the Fourier tranform of input quantities. Options are “fft”, for a standard Fast Fourier Transform as implemented by numpy.fft, and “sinc”, for the sinc function decomposition described in Uttley and Malzac 2023 (arXiv:2312.08302).

Attributes:

times: np.array(float)

The array of times over which quantities are to be Fourier transformed. Necessary to utilize the sinc transform method. This array needs to be linearly spaced to use the fft method; if this is not the case, the class defaults to the sinc method.

time_bins: np.array(float)

The array of widths of each bin in the “time” array. The FFT method can only be called if all the elements of the array are equal to the first, or in other words if the “time” array is linearly spaced.

n_times: int

The size of the “times” and “time_bins” arrays.

method: string{“fft” | “sinc” }

The computational method to calculate the Fourier tranform of input quantities. Options are “fft”, for a standard Fast Fourier Transform as implemented by numpy.fft, and “sinc”, for the sinc function decomposition described in Uttley and Malzac 2023 (arXiv:2312.08302).

freqs: np.array(float)

The array of Fourier frequencies over which the Fourier transform is computed. When using the “fft” method, this is derived from the “times” array, considering only the positive frequency bins. When using the “sinc” method, this is set by the user when the operator object is first initialized.

n_freqs: int

The size of the “freqs” array.

irf_sinc_arr: np.array(complex,complex)

A matrix of size (n_freqs) times n_times, which stores how each time bin in the “times” array maps to a given sinc function with frequency given in the “freqs” array. This is required for use of the sinc method, and therefore is updated automatically every time the user changes frequency grid.

transform(input_array)

This method acts like a wrapper to Fourier transform array, and calls the appropriate method depending on user settings.

Parameters:

input_array: np.array(float)

The array to be Fourier transformed.

Returns:

transform: np.array(complex)

The Fourier transform of the input array.

class ndspec.Timing.PowerSpectrum(times, freqs=0, method='fft', verbose=False)

This class computes model power spectra from a given signal input. The public methods in the class call those in nDspecOperator and FourierProduct nas ecessary to Fourier transform input quantities and handle changes to the internal frequency grid.

Parameters inherited from FourierProduct:

times: np.array(float)

The array of times over which the input signal is defined.

freqs: np.array(float)

The Fourier frequency array over which the power spectrum is defined.

method: string{“fft” | “sinc” }, default=’fft’

The computational method to calculate the power spectrum. Options are “fft”, for a standard Fast Fourier Transform as implemented by numpy.fft, and “sinc”, for the sinc function decomposition described in Uttley and Malzac 2023 (arXiv:2312.08302).

verbose: bool

A setting to display initialization warnings.

Attributes inherited from FourierProduct:

times: np.array

The array of times over which the input signal is defined.

time_bins: np.array

The array of widths of each bin in the “time” array. The FFT method can only be used if all the elements of the array are equal to the first, or in other words if the “time” array is linearly spaced.

n_times: int

The size of the “times” and “time_bins” arrays.

method: {“fft” | “sinc” }

The computational method to calculate the Fourier tranform of input quantities. Options are “fft”, for a standard Fast Fourier Transform as implemented by numpy.fft, and “sinc”, for the sinc function decomposition described in Uttley and Malzac 2023 (arXiv:2312.08302).

freqs: np.array(float)

The array of Fourier frequencies over which the power spectrum is defined and/or computed.

n_freqs: int

The size of the “freqs” array.

irf_sinc_arr: np.array(complex,complex)

A matrix of size (n_freqs x n_times), which stores how each time bin in the “times” array maps to a given sinc function with frequency given in the “freqs” array. This is required for use of the sinc method, and therefore is updated automatically every time the user changes frequency grid.

Other Attributes:

power_spec: np.array(float)

An array of size(n_freqs), storing the un-normalized power spectrum of an input signal.

compute_psd(signal=None, freq=None, params=None)

This method calculates the power spectrum, defined over the class “freqs” Fourier frequency array, of an input array, and assigns it to the power_spec attribute contained in the class instance.

Parameters:

signal: np.array(float), default None

The quantity from which to calculate the power spectrum. If it is not provided, the class will use the model object stored in the ``model’’ attribute instead.

freq: np.array(float), default None

The array of Fourier frequencies over which to compute the model. If none is passed and evaluating a Fourier-domain model, uses the frequency stored in the object at initialization.

params: lmfit.Parameters, default None

The parameter values to use in evaluating the model. If none are provided, the ‘’model_params’’ attribute stored in the object is used.

plot_psd(units='Power*freq', return_plot=False)

This method plots the either the power or power per unit frequency as a function of frequency, stored in the class instance.

Parameters:

units: string , default=”Power*freq”

Sets the units to plot on the y axis - the default “power*freq” displays the power at that frequency, “power” instead uses the power per unit frequency.

Returns:

fig: matplotlib.figure, optional

The plot object produced by the method.

rebin_frequency(new_grid, use_log=True)

This methods rebins the class freqs array, and also updates all relevant attributes (n_freqs, power_spec). The re-binned power spectrum is calculated by interpolating the previous power over the new grid.

Parameters:

new_grid: np.array(float)

The new grid of Fourier frequencies over which to rebin the power spectrum.

Other parameters:

use_log: bool, default True

Switches between interpolating the power spectrum, or the base 10 logarithm of the power spectrum, which is more accurate if the power varies by many orders of magnitude with frequency

set_psd_model(model, params=None)

This method sets a LMFit CompositeModel object to be used to evaluate the power spectrum; useful if you want your model to be defined in the Fourier rather than time domain.

Parameters:

model: lmfit.model or lmfit.compositemodel

The lmfit wrapper of the model one wants to fit to the data.

params: lmfit.Parameters, default: None

The parameter values from which to start evalauting the model during the fit. If it is not provided, all model parameters will default to 0, set to be free, and have no minimum or maximum bound.

class ndspec.Timing.CrossSpectrum(times, freqs=0, energ=None, method='fft', verbose=False)

This class computes model cross spectra from a given input, assumed to be either a user-provided impulse response function or a user-provided model (such as a transfer function) already defined in Fourier space. It can handle both one-dimensional, frequency dependent cross spectra between just two energy bands, or two-dimensional cross spectra defined over multiple energy and frequency grids.

Parameters inherited from FourierProduct:

times: np.array(float)

The array of times over which the input signal is defined.

freqs: np.array(float)

The Fourier frequency array over which the power spectrum is defined.

method: string{“fft” | “sinc” }, default=”fft”

The computational method to calculate the cross spectrum. Options are “fft”, for a standard Fast Fourier Transform as implemented by numpy.fft, and “sinc”, for the sinc function decomposition described in Uttley and Malzac 2023 (arXiv:2312.08302).

verbose: bool

A setting to trigger initialization warnings.

Other parameters:

energ: np.array(float)

The array of energy (or energy channels) bin mid-points over which a two-dimensional cross spectrum is defined. If we are only considering the cross-spectrum between two energy bands, this defaults to “None”.

Attributes inherited from FourierProduct:

times: np.array(float)

The array of times over which the input signal is defined.

time_bins: np.array(float)

The array of widths of each bin in the “time” array. The FFT method can only be used if all the elements of the array are equal to the first, or in other words if the “time” array is linearly spaced.

n_times: int

The size of the “times” and “time_bins” arrays.

method: string{“fft” | “sinc” }

The computational method to calculate the Fourier tranform of input quantities. Options are “fft”, for a standard Fast Fourier Transform as implemented by numpy.fft, and “sinc”, for the sinc function decomposition described in Uttley and Malzac 2023 (arXiv:2312.08302).

freqs: np.array(float)

The array of Fourier frequencies over which the power spectrum is defined and/or computed.

n_freqs: int

The size of the “freqs” array.

irf_sinc_arr: np.array(complex,complex)

A matrix of size (n_freqs x n_times), which stores how each time bin in the “times” array maps to a given sinc function with frequency given in the “freqs” array. This is required for use of the sinc method, and therefore is updated automatically every time the user changes frequency grid.

Other attributes:

n_chans: int

The size of the “energy” array. Defaults to 1 for a one-dimensional cross spectrum between two energy bands.

chans: np.array(int)

An array containing the indexes of each energy bin (or channel) in the cross spectrum. Defaults to 0 for a one-dimensional cross spectrum between two energy bands.

power_spec: np.array(float)

A one-dimensional array of size (n_freqs) containing the assumed shape of the power spectrum. This is necessary to a) calculate the cross product from an impulse response and b) calculate Fourier frequency averaged products like lag-energy spectra in a given frequency range.

imp_resp: np.array(float,float)

An array of size (n_chans x n_times), storing the model impulse response function from which to compute the cross spectrum.

ref: np.array(float)

An array of size (n_times), storing the reference band signal. The reference band can either be computed directly from the imp_resp attribute, using the set_reference_idx method, or it can be provided by the user, using the set_reference_lc method.

correct_ref: bool

A bool that sets whether the reference band used in the cross spectrum calculations is corrected by subtracting the channel of interest (correct_ref=True) or if it is kept the same (correct_ref=False).

trans_func: np.array(complex,complex)

An array of size (n_chans x n_freqs), containing the Fourier transform of each channel of the impulse response function provided by the user - formally, this is the transfer function. It is used together with the weighing power spectrum power_spec to calculate the cross spectrum, and is necessary to compute one-dimensional frequency-dependent products (such as lag-frequency spectra) from a two-dimensional cross spectrum. If users are not providing an impulse response function, this attribute can be set by hand (e.g. from a model defined in Fourier space).

cross: np.array(complex,complex)

An array of size (n_chans x n_freqs) containing the full (one or two dimensional) cross spectrum, defined as the cross product between the Fourier transforms of one or more channels of interest and of the reference band, and weighed by the power spectrum. If users are not providing an impulse response or transfer function, this attribute can be set by hand (e.g. from a model defined in Fourier space).

cross_from_irf(signal=None, ref_bounds=None, power=None)

This method computes the transfer function and cross spectrum from the impulse response, reference band, and power spectra provided by the user. These can either be already stored by the setter methods set_psd_weights, set_impulse, and set_reference_idx or set_reference_lc (which is the default behavior), or they can be passed as arguments of this method. In the latter case, the reference can only be provided in array, rather than channel index, form.

Parameters:

signal: np.array(float,float), default=self.imp_resp

An array of size (n_chans x n_times) containing the model impulse response function.

ref_bounds: np.array(float)

A list with lower and upper energy channel bounds to be used in the reference band. By default, we assume that the reference band is identical to that used in calculating the cross spectrum. As implemented here, the limits specified in ref_bounds are included in the reference - e.g. [0.3,10.0] keV rather than (0.3,10.) keV.

power: np.array(float) or PowerSpectrum, default=self.power_spec

Either an array of size (n_freqs), or a PowerSpectrum nDspec object, that is to be used as the weighing power spectrum when computing the cross spectrum.

cross_from_transfer(signal=None, ref_bounds=None, power=None)

This method computes thecross spectrum from an input defined in Fourier space (such as a transfer function), reference band, and power spectra provided by the user. These can either be already stored by the setter methods set_psd_weights, set_transfer, and set_reference_idx or set_reference_lc (which is the default behavior), or they can be passed as arguments of this method. In the latter case, the reference can only be provided in array, rather than channel index, form.

Parameters:

signal: np.array(float,float), default=self.trans_func

An array of size (n_chans x n_freqs) containing a model transfer function defined in Fourier space.

ref_bounds: np.array(float)

A list with lower and upper energy channel bounds to be used in the reference band. By default, we assume that the reference band is identical to that used in calculating the cross spectrum. As implemented here, the limits specified in ref_bounds are included in the reference - e.g. [0.3,10.0] keV rather than (0.3,10.) keV.

power: np.array(float) or PowerSpectrum, default=self.power_spec

Either an array of size (n_freqs), or a PowerSpectrum nDspec object, that is to be used as the weighing power spectrum when computing the cross spectrum.

imag()

This method returns the imaginary part of the cross spectrum, both for one and two dimensional cases.

Returns:

imag: np.array(float)

An array of imaginary parts of the cross spectrum, of size (n_chans x n_freqs)

imag_energy(freq_bounds)

This method computes the imaginary part of a one-dimensional cross spectrum as a function of energy, by averaging the attribute “cross” over a given frequency range specified by the user.

Parameters:

freq_bounds: np.array(float)

A list with lower and upper frequency bounds over which to average the two dimensinoal cross spectrum

Returns:

imag_spectrum: np.array(float)

An array of size (n_chans) containing the imaginary part of the Fourier frequency averaged cross spectrum, as a function of energy.

imag_frequency(int_bounds, ref_bounds=None)

This method computes the imaginary part of a one-dimensional, Fourier frequency dependent cross spectrum, by summing over energy channels specified by the user and (if a new reference band is used) re-calculating the cross product between the reference band and channels of interest.

Parameters:

int_bounds: np.array(float)

A list with lower and upper energy channel bounds to be used in the channels of interest - using the convention used in X-ray spectral timing, this should be the energy band where reverberation lags appear with negative values.

ref_bounds: np.array(float), default=None

A list with lower and upper energy channel bounds to be used in the reference band. By default, we assume that the reference band is identical to that used in calculating the cross spectrum.

Returns:

imag_spectrum: np.array(float)

A one-dimensional array of size (n_freqs) containing the imaginary part of the one dimensional cross spectrum between the channels of interest and the reference band provided by the user, as a function of Fourier frequency.

lag()

This method converts the phase of the cross spectrum, both for one and and two dimensional cases, into time lags.

Returns:

lags: np.array(float)

An array of phases of the cross spectrum, of size (n_chans x n_freqs)

lag_energy(freq_bounds)

This method computes the time lags of a one-dimensional cross spectrum as a function of energy, by averaging the attribute “cross” over a given frequency range specified by the user.

Parameters:

freq_bounds: np.array(float)

A list with lower and upper frequency bounds over which to average the two dimensinoal cross spectrum

Returns:

lag_spectrum: np.array(float)

An array of size (n_chans) containing the time lags of the Fourier frequency averaged cross spectrum, as a function of energy.

lag_frequency(int_bounds, ref_bounds=None)

This method computes the time lag of a one-dimensional, Fourier frequency dependent cross spectrum, by summing over energy channels specified by the user and (if a new reference band is used) re-calculating the cross product between the reference band and channels of interest.

Parameters:

int_bounds: np.array(float)

A list with lower and upper energy channel bounds to be used in the channels of interest - using the convention used in X-ray spectral timing, this should be the energy band where reverberation lags appear with negative values.

ref_bounds: np.array(float), default=None

A list with lower and upper energy channel bounds to be used in the reference band. By default, we assume that the reference band is identical to that used in calculating the cross spectrum.

Returns:

mod_spectrum: np.array(float)

A one-dimensional array of size (n_freqs) containing the time lags of the one dimensional cross spectrum between the channels of interest and the reference band provided by the user, as a function of Fourier frequency.

mod()

This method returns the modulus of the cross spectrum, both for one and and two dimensional cases.

Returns:

mod: np.array(float)

An array of moduli of the cross spectrum, of size (n_chans x n_freqs)

mod_energy(freq_bounds)

This method computes the modulus of a one-dimensional cross spectrum as a function of energy, by averaging the attribute “cross” over a given frequency range specified by the user.

Parameters:

freq_bounds: np.array(float)

A list with lower and upper frequency bounds over which to average the two dimensinoal cross spectrum

Returns:

mod_spectrum: np.array(float)

An array of size (n_chans) containing the modulus of the Fourier frequency averaged cross spectrum, as a function of energy.

mod_frequency(int_bounds, ref_bounds=None)

This method computes the modulus of a one-dimensional, Fourier frequency dependent cross spectrum, by summing over energy channels specified by the user and (if a new reference band is used) re-calculating the cross product between the reference band and channels of interest.

Parameters:

int_bounds: np.array(float)

A list with lower and upper energy channel bounds to be used in the channels of interest - using the convention used in X-ray spectral timing, this should be the energy band where reverberation lags appear with negative values.

ref_bounds: np.array(float), default=None

A list with lower and upper energy channel bounds to be used in the reference band. By default, we assume that the reference band is identical to that used in calculating the cross spectrum.

Returns;

mod_spectrum: np.array(float)

A one-dimensional array of size (n_freqs) containing the modulus of the one dimensional cross spectrum between the channels of interest and the reference band provided by the user.

phase()

This method returns the phase of the cross spectrum, both for one and and two dimensional cases.

Returns:

phases: np.array(float)

An array of phases of the cross spectrum, of size (n_chans x n_freqs)

phase_energy(freq_bounds)

This method computes the phase of a one-dimensional cross spectrum as a function of energy, by averaging the attribute “cross” over a given frequency range specified by the user.

Parameters:

freq_bounds: np.array(float)

A list with lower and upper frequency bounds over which to average the two dimensinoal cross spectrum

Returns:

phase_spectrum: np.array(float)

An array of size (n_chans) containing the phase of the Fourier frequency averaged cross spectrum, as a function of energy.

phase_frequency(int_bounds, ref_bounds=None)

This method computes the phase of a one-dimensional, Fourier frequency dependent cross spectrum, by summing over energy channels specified by the user and (if a new reference band is used) re-calculating the cross product between the reference band and channels of interest.

Parameters:

int_bounds: np.array(float)

A list with lower and upper energy channel bounds to be used in the channels of interest - using the convention used in X-ray spectral timing, this should be the energy band where reverberation lags appear with negative values.

ref_bounds: np.array(float), default=None

A list with lower and upper energy channel bounds to be used in the reference band. By default, we assume that the reference band is identical to that used in calculating the cross spectrum.

Returns:

phase_spectrum: np.array(float)

A one-dimensional array of size (n_freqs) containing the phase one dimensional cross spectrum between the channels of interest and the reference band provided by the user, as a function of Fourier frequency.

plot_cross_1d(form='polar', return_plot=False)

This method plots the a one-dimensional cross spectrum as a function of Fourier frequency.

Parameters:

form: string, default=”polar”

A qualifier to choose in which units to plot the cross spectrum. By default, form=”polar” will plot the modulus, phase, and lag frequency spectrum. Alternatively, form=”cartesian” plots the real and imaginary parts of the cross spectrum.

Returns:

fig: matplotlib.figure, optional

The plot object produced by the method.

plot_cross_2d(form='polar', energy_limits=[0.3, 10.5], return_plot=False, normalize_en=True)

Plots the a two-dimensional cross spectrum as a function of Fourier frequency and energy.

Parameters:

form: string,default=”polar”

A qualifier to choose in which units to plot the cross spectrum. By default, form=”polar” will plot the modulus, phase, and lag spectrum. Alternatively, form=”cartesian” plots the real and imaginary parts of the cross spectrum.

energy_limits: list(float)

The lower and upper bound to be used in the energy axis of the cross spectrum.

Returns:

fig: matplotlib.figure, optional

The plot object produced by the method.

real()

This method returns the real part of the cross spectrum, both for one and two dimensional cases.

Returns:

real: np.array(float)

An array of real parts of the cross spectrum, of size (n_chans x n_freqs)

real_energy(freq_bounds)

This method computes the real part of a one-dimensional cross spectrum as a function of energy, by averaging the attribute “cross” over a given frequency range specified by the user.

Parameters:

freq_bounds: np.array(float)

A list with lower and upper frequency bounds over which to average the two dimensinoal cross spectrum

Returns:

real_sectrum: np.array(float)

An array of size (n_chans) containing the real part of the Fourier frequency averaged cross spectrum, as a function of energy.

real_frequency(int_bounds, ref_bounds=None)

This method computes the real part of a one-dimensional, Fourier frequency dependent cross spectrum, by summing over energy channels specified by the user and (if a new reference band is used) re-calculating the cross product between the reference band and channels of interest.

Parameters:

int_bounds: np.array(float)

A list with lower and upper energy channel bounds to be used in the channels of interest - using the convention used in X-ray spectral timing, this should be the energy band where reverberation lags appear with negative values.

ref_bounds: np.array(float), default=None

A list with lower and upper energy channel bounds to be used in the reference band. By default, we assume that the reference band is identical to that used in calculating the cross spectrum.

Returns:

real_spectrum: np.array(float)

A one-dimensional array of size (n_freqs) containing the real part of he one dimensional cross spectrum between the channels of interest and the reference band provided by the user, as a function of Fourier frequency.

rebin_frequency(new_grid)

This method rebins the cross spectrum and transfer function stored in the object to a new frequency grid, and updates the relevant class attributes as necessary.

Parameters:

new_grid: np.array(float)

An array of arbitrary size, containing the Fourier frequency bin midpoints of the new grid.

set_impulse(signal)

This method sets the impulse response function imp_resp from which the cross spectrum can be calculated.

Parameters:

signal: np.array(float,float)

An array of size (n_chans x n_times) containing the model impulse response function.

set_psd_weights(input_power)

This method sets the weighing power spectrum power_spec from a given input.

Parameters:

input_power: np.array(float) or PowerSpectrum

Either an array of size (n_freqs) that is to be used as the weighing power spectrum when computing the cross spectrum, or an nDspec PowerSpectrum object. Both have to be defined over the same Fourier frequency array.

set_reference_energ(ref_bounds, correct_ref=True)

This method sets the reference band from a range of energies provided by the user. Users can also specify whether or not they want each channel of interest to be subtracted from the reference band when computing the cross spectrum.

Parameters:

ref_bounds: np.array(float)

A list with lower and upper energy channel bounds to be used in the reference band. By default, we assume that the reference band is identical to that used in calculating the cross spectrum. As implemented here, the limits specified in ref_bounds are included in the reference - e.g. [0.3,10.0] keV rather than (0.3,10.) keV.

correct_ref: bool, default=True

Flag to correct the reference band by removing each channel of interest or not.

set_reference_lc(input_lc, correct_ref=False)

This method sets the reference band from an array (e.g. count rate as a function of energy) provided by the user. Users can also specify whether or not they want each channel of interest to be subtracted from the reference band when computing the cross spectrum.

Parameters:

input_lc: np.array(float)

An arrray of size (n_chans) containing either the reference band lightcurve (if the model is defined in the time domain) or its Fourier transform (for models defined in the Fourier domain).

correct_ref: bool, default=False

Flag to correct the reference band by removing each channel of interest or not.

set_transfer(signal)

This method sets the transfer function trans_func from which the cross spectrum can be calculated.

Parameters:

signal: np.array(float,float)

An array of size (n_chans x n_freqs) containing the model transfer function.

transfer_from_irf(signal=None)

This method computes and store the transfer function from a given impulse response provided by the user. This can either be already stored by the setter method set_impulse (which is the default behavior), or it can be passed as argument of this method.

Parameters:

signal: np.array(float,float), default=self.imp_resp

An array of size (n_chans x n_times) containing the model impulse response function.

SimpleFit Classes

class ndspec.SimpleFit.SimpleFit(likelihood='chisq')

Generic least-chi squared fitter class, used internally to store methods that are shared between all the fitter types.

Attributes:

model: lmfit.CompositeModel

A lmfit CompositeModel object, which contains a wrapper to the model component(s) one wants to fit to the data.

model_params: lmfit.Parameters

A lmfit Parameters object, which contains the parameters for the model components.

likelihood: str

A string that allows to switch between different fit statistics; which one is available depends on the type of fitter object. Uses chi-squared likelihood by default. Users can set different likelihoods with the appropriate setter method.

custom_likelihood: function

A function users can set to bypass the supported likelihoods and instead provide their own.

custom_args: tuple

A tuple including any custom arguments (in addition to the data and model values to be compared) necessary to calculate the custom likelihood

fit_result: lmfit.MinimizeResult

A lmfit MinimizeResult, which stores the result (including best-fitting parameter values, fit statistics etc) of a fit after it has been run.

data: np.array(float)

An array storing the data to be fitted. If the data is complex and/or multi-dimensional, it is flattened to a single dimension in order to be compatible with the LMFit fitter methods.

data_err: np.array(float)

An array containing the uncertainty on the data to be fitted. It is also stored as a one-dimensional array regardless of the type or dimensionality of the initial data.

noise: np.array(float) or None

If loaded, an array containing the background spectrum, including only the channels noticed in the fit.

noise_err: np.array(float or None)

If loaded, an array containing the sqrt of the background counts, only in the channels noticed during the fit. Used to compute the fit statistic.

_data_unmasked, _data_err_unmasked, _noise_unmasked: np.array(float)

The arrays of every data bin, its error and (if loaded) the backgruond, regardless of which ones are ignored or noticed during the fit. Used exclusively to enable book keeping internal to the fitter class.

fit_data(algorithm='leastsq')

This method attempts to minimize the residuals of the model with respect to the data defined by the user. The fit always starts from the set of parameters defined with .set_params(). Once the algorithm has completed its run, it prints to terminal the best-fitting parameters, fit statistics, and simple selection criteria (reduced chi-squared, Akaike information criterion, and Bayesian information criterion).

Parameters:

algorithm: str, default=”leastsq”

The fitting algorithm to be used in the minimization. The possible choices are detailed on the LMFit documentation page: https://lmfit.github.io/lmfit-py/fitting.html#fit-methods-table.

get_residuals(res_type, model=None, mask=True)

This methods return the residuals (either as data/model, or as contribution to the total chi squared) of the input model, given the parameters set in model_parameters, with respect to the data.

Parameters:

res_type: string

If set to “ratio”, the method returns the residuals defined as data/model. If set to “chisq”, it returns the contribution of each energy channel to the total chi squared. If set to “custom”, the residuals are based on whatever custom likelihood the user defined.

mask: bool, default True

A flag to decide whether to compare the model against the masked or unmasked data.

Returns:

residuals: np.array(float)

An array of the same size as the data, containing the model residuals in each channel.

bars: np.array(float)

An array of the same size as the residuals, containing the one sigma range for each contribution to the residuals.

print_fit_report()

This method prints the current fit result.

print_fit_stat()

This method compares the model defined by the user, using the last set of parameters to have been set in the class, to the data stored. It then prints the chi-squared goodness-of-fit to terminal, along with the number of data bins, free parameters and degrees of freedom.

print_model()

This method prints out model components, model parameters, and their settings.

set_custom_likelihood(likelihood_function, *args)

This method allows users to define their own custom likelihood function, which can then optimized during a fit. In addition, this sets the value of the class “likelihood” string to custom, to signal to the other methods that a custom likelihood is in use and should be used for plots, residuals etc.

Parameters:

likelihood_function: function

The name of the function which calculates the model residuals; e.g., if we want to minimize the difference between data and model, we would define: def diff(data,model):

return data-model

and call set_custom_likelihood(diff).

*args:

Additional arguments to be passed to the likelihood calculation, excluding the data and model (which are always included automatically by the class). Following the example above: def diff(data,model,factor):

return factor*(data-model)

and call set_custom_likelihood(diff,5) - if we want to set “factor” to 5.

set_model(model, params=None)

This method is used to pass the model users want to fit to the data. Optionally it is also possible to pass the initial parameter values of the model.

Parameters:

model: lmfit.model or lmfit.compositemodel

The lmfit wrapper of the model one wants to fit to the data.

params: lmfit.Parameters, default: None

The parameter values from which to start evalauting the model during the fit. If it is not provided, all model parameters will default to 0, set to be free, and have no minimum or maximum bound.

set_params(params)

This method is used to set the model parameter names and values. It can be used both to initialize a fit, and to test different parameter values before actually running the minimization algorithm.

Parameters:

params: lmfit.parameter

The parameter values from which to start evalauting the model during the fit.

class ndspec.SimpleFit.EnergyDependentFit

Internal book-keeping class used to manage noticing or ignoring energy channels, for cases when the data requires an instrument response.

Stores the full (unmasked) energy center/bounds, and data arrays, a mask used to track which channels/data points are noticed or ignored, as well as the masked arrays containing only the noticed bins.

Attributes:

energs: np.array(float)

The array of physical photon energies over which the model is computed. Defined as the middle of each bin in the energy range stored in the instrument response provided.

energ_bounds: np.array(float)

The array of energy bin widths, for each bin over which the model is computed. Defined as the difference between the uppoer and lower bounds of the energy bins stored in the insrument response provided.

ear: np.array(float)

The array of energy bin bounds, for each bin over which the model is computed. Only necessary when calling Xspec models due to their unique input structure.

ebounds: np.array(float)

The array of energy channel bin centers for the instrument energy channels, as stored in the instrument response provided. Only contains the channels that are noticed during the fit.

ewidths: np.array(float)

The array of energy channel bin widths for the instrument energy channels, as stored in the instrument response provided. Only contains the channels that are noticed during the fit.

ebounds_mask: np.array(bool)

The array of instrument energy channels that are either ignored or noticed during the fit. A given channel i is noticed if ebounds_mask[i] is True, and ignored if it is false.

n_chans: int

The number of channels that are to be noticed during the fit.

_all_chans: int

The total number of channels in the loaded response matrix.

n_bins: int

Only used for two-dimensional data fitting. Defined as the number of noticed channels, times the number of bins in the second dimension (e.g. Fourier frequency).

_all_bins: int

Only used for two-dimensional data fitting. Defined as the total number of channels, times the number of bins in the second dimension (e.g. Fourier frequency).

_emin_unmasked, _emax_unmasked, _ebounds_unmasked, _ewidths_unmasked: np.array(float)

The array of every lower bound, upper bound, channel center and channel widths stored in the response, regardless of which ones are ignored or noticed during the fit. Used exclusively to facilitate book-keeping internal to the fitter class.

ignore_energies(bound_lo, bound_hi)

This method adjusts the arrays stored such that they (and the fit) ignore selected channels based on their energy bounds.

Parameters:

bound_lofloat

Lower bound of ignored energy interval.

bound_hifloat

Higher bound of ignored energy interval.

notice_energies(bound_lo, bound_hi)

This method adjusts the data arrays stored such that they (and the fit) notice selected (previously ignore) channels based on their energy bounds.

Parameters:

bound_lofloat

Lower bound of ignored energy interval.

bound_hifloat,

Higher bound of ignored energy interval.

class ndspec.SimpleFit.FrequencyDependentFit(freqs)

Internal book-keeping class used to manage noticing or ignoring Fourier frequency bins.

Stores the full (unmasked) Fourier bins, and data arrays, a mask used to track which bins/data points are noticed or ignored.

Attributes:

_freqs_unmasked: np.array(float)

If the data and model explicitely depend on Fourier frequency (e.g. a power spectrum), this is the array of Fourier frequency over which all data and model are defined, including bins that are ignored in the fit.

If instead the data depends from some other energy (e.g. energy), it contains both noticed and ignored frequency intervals over which to produce spectral-timing products. For example, a user might input a set of 7 ranges of frequencies to calculate lag energy spectra, but only want to consider the first and last 3, and ignore the middle one.

freqs_mask np.array(bool)

The array of Fourier frequencies that are either ignored or noticed during the fit. A given channel i is noticed if freqs_mask[i] is True, and ignored if it is false.

n_freqs: int

The number of Fourier frequency bins that are noticed in the fit.

n_bins: int

Only used for two-dimensional data fitting. Defined as the number of noticed channels, times the number of bins in the second dimension (e.g. Fourier frequency).

_all_bins: int

Only used for two-dimensional data fitting. Defined as the total number of channels, times the number of bins in the second dimension (e.g. Fourier frequency).

ignore_frequencies(bound_lo, bound_hi)

This method adjusts the arrays stored such that they (and the fit) ignore selected frequencies based on user-supplied bounds bounds.

Parameters:

bound_lofloat

Lower bound of ignored frequency interval.

bound_hifloat

Higher bound of ignored frequency interval.

notice_frequencies(bound_lo, bound_hi)

This method adjusts the arrays stored such that they (and the fit) ignore selected frequencies based on user-supplied bounds bounds.

Parameters:

bound_lofloat

Lower bound of ignored frequency interval.

bound_hifloat

Higher bound of ignored frequency interval.

Data loading utilities

ndspec.SimpleFit.load_lc(path)

This function loads an X-ray lightcurve, given an input path to an OGIP-compatible file.

Parameters:

path: str

A string pointing to the lightcurve file to be loaded

Returns:

time_bins: np.array(float)

An array of time stamps covered by the lightcurve

counts: np.array(float)

An array of counts rates (defined in counts per second) contained in the lightcurve

gti: list([float,float])

A list of good time intervals over which the lightcurve is defined.

ndspec.SimpleFit.load_pha(path, response)

This function loads an X-ray spectrum , given an input path to an OGIP-compatible file and a nDspec ResponseMatrix object to be applied to the spectrum.

Parameters:

path: str

A string pointing to the spectrum file to be loaded

response: nDspec.ResponseMatrix

The instrument response matrix, loaded in nDspec, corresponding to the spectrum to be loaded

Returns:

bin_bounds_lo: np.array(float)

An array of lower energy channel bounds, in keV, as contained in the input file. If the spectrum was grouped, this contains the lower bounds of the spectrum after rebinning.

bin_bounds_hi: np.array(float)

An array of upper energy channel bounds, in keV, as contained in the input file. If the spectrum was grouped, this contains the lower bounds of the spectrum after rebinning.

counts_per_group: np.array(float)

The total number of photon counts in each energy channel. If the spectrum was grouped, this contains the counts in each channel after rebinning.

spectrum_error: np.array(float)

The error on the counts in each group, including both Poisson and (if present) systematic errors

exposure: float

The exposure time contained in the spectrum file.

backscal: float

The background scaling factor. Typically used by imaging instruments to account for different extraction region size for the source and background.

FitPowerSpectrum Class

class ndspec.FitPowerSpectrum.FitPowerSpectrum(likelihood='chisq')

Least-chi squared fitter class for a powerspectrum, defined as the product between a Fourier-transformed lightcurve and its complex conjugate.

Given an array of Fourier frequencies, a power spectrum, its error and a model (defined in Fourier space), this class handles fitting internally using the lmfit library.

Poisson noise in the data is not accounted for explicitely. Users should either pass noise-subtracted data, or include a constant component in the model definition (see below).

Attributes inherited from SimpleFit:

model: lmfit.CompositeModel

A lmfit CompositeModel object, which contains a wrapper to the model component(s) one wants to fit to the data.

model_params: lmfit.Parameters

A lmfit Parameters object, which contains the parameters for the model components.

likelihood: str

A string that allows to switch between different fit statistics; which one is available depends on the type of fitter object. Uses chi-squared likelihood by default. Users can set different likelihoods either at initialization or with the appropriate setter method.

custom_likelihood: function

A function users can set to bypass the supported likelihoods and instead provide their own.

custom_args: tuple

A tuple including any custom arguments (in addition to the data and model values to be compared) necessary to calculate the custom likelihood

fit_result: lmfit.MinimizeResult

A lmfit MinimizeResult, which stores the result (including best-fitting parameter values, fit statistics etc) of a fit after it has been run.

data: np.array(float)

An array storing the data to be fitted. Only contains noticed bins.

data_err: np.array(float)

An array containing the uncertainty on the data to be fitted. Only contains noticed bins.

_data_unmasked, _data_err_unmasked: np.array(float)

The arrays of every data bin and its error, regardless of which ones are ignored or noticed during the fit. Used exclusively to enable book keeping internal to the fitter class.

Attributes inherited from FrequencyDependentFit:

_freqs_unmasked

freqs_mask

n_freqs

Other attributes:

freqs: np.array(float)

The Fourier frequency over which both the data and model are defined, in units of Hz. Only contains noticed bins.

eval_model(params=None, freq=None, mask=True)

This method is used to evaluate and return the model values for a given set of parameters, over a given Fourier frequency grid. By default it will evaluate the model over the data Fourier frequency grid and using the values stored internally in the model_params attribute, but passing a different grid/set of parameters is also possible.

Parameters:

params: lmfit.Parameters, default None

The parameter values to use in evaluating the model. If none are provided, the model_params attribute is used.

freq: np.array(float), default None

The the Fourier frequencies over which to evaluted the model. If none are provided, the same frequencies over which the data is defined are used.

mask: bool, default True

A boolean switch to choose whether to mask the model output to only include the noticed energy channels, or to also return the ones that have been ignored by the users.

Returns:

model: np.array(float)

The model evaluated over the given Fourier frequency array, for the given input parameters.

plot_data(units='fpower', return_plot=False)

This method plots the powerspectrum loaded by the user as a function of Fourier frequency. It is possible to plot both units of power, and power times frequency, depending on user input.

It is also possible to return the figure object, for instance in order to save it to file.

Parameters:

units: str, default=”fpower”

The units to use for the y axis. units=”fpower”, the default, plots the data in units of power*frequency. units=”power” instead plots the data in units of power.

return_plot: bool, default=False

A boolean to decide whether to return the figure objected containing the plot or not.

Returns:

fig: matplotlib.figure, optional

The plot object produced by the method.

plot_model(plot_data=True, plot_components=False, params=None, units='fpower', residuals='chisq', return_plot=False)

This method plots the model defined by the user as a function of Fourier frequency, as well as (optionally) its components, and the data plus model residuals. It is possible to plot both units of power, and power times frequency, depending on user input.

It is also possible to return the figure object, for instance in order to save it to file.

Parameters:

plot_data: bool, default=True

If true, both model and data are plotted; if false, just the model.

plot_components: bool, default=False

If true, the model components are overplotted; if false, they are not. Only additive model components will display their values correctly.

params: lmfit.parameters, default=None

The parameters to be used to evaluate the model. If False, the set of parameters stored in the class is used

units: str, default=”fpower”

The units to use for the y axis. units=”fpower”, the default, plots the data in units of power*frequency. units=”power” instead plots the data in units of power.

residuals: str, default=”chisq”

The units to use for the residuals. If residuals=”chisq”, the plot shows the residuals in units of data-model/error; if residuals=”ratio”, the plot instead uses units of data/model. If a custom likelihood is defined, then the residuals are computed from that likelihood.

return_plot: bool, default=False

A boolean to decide whether to return the figure objected containing the plot or not.

Returns:

fig: matplotlib.figure, optional

The plot object produced by the method.

set_data(data, data_err=None, data_grid=None)

This method is used to pass the data users want to fit. The input can be either three arrays including the power, its erorr, and the Fourier frequency grid, or a Stingray Powerspectrum object, in which case the error and frequency array are handled automatically.

Parameters:

data: np.array(float) or stingray.powerspectrum

The power spectrum to be fitted, either as a numpy array or a stingray object

data_err: np.array(float), default: None

The error on the power spectrum. If passing a stingray object, this is not necessary and is therefore ignored

data_grid: np.array(float), default: None

The Fourier frequency grid over which the data (and model) are defined. If passing a stingray object, this is not necessary and is therefore ignored

FitTimeAvgSpectrum Class

class ndspec.FitTimeAvgSpectrum.FitTimeAvgSpectrum(likelihood='chisq')

Least-chi squared fitter class for a time averaged spectrum, defined as the count rate as a function of photon channel energy bound.

Given an instrument response, a count rate spectrum, its error and a model (defined in energy space), this class handles fitting internally using the lmfit library.

Attributes inherited from SimpleFit:

model: lmfit.CompositeModel

A lmfit CompositeModel object, which contains a wrapper to the model component(s) one wants to fit to the data.

model_params: lmfit.Parameters

A lmfit Parameters object, which contains the parameters for the model components.

likelihood: str

A string that allows to switch between different fit statistics; which one is available depends on the type of fitter object. Uses chi-squared likelihood by default. Users can set different likelihoods either at initialization or with the appropriate setter method.

custom_likelihood: function

A function users can set to bypass the supported likelihoods and instead provide their own.

custom_args: tuple

A tuple including any custom arguments (in addition to the data and model values to be compared) necessary to calculate the custom likelihood

fit_result: lmfit.MinimizeResult

A lmfit MinimizeResult, which stores the result (including best-fitting parameter values, fit statistics etc) of a fit after it has been run.

data: np.array(float)

An array storing the data to be fitted. If the data is complex and/or multi-dimensional, it is flattened to a single dimension in order to be compatible with the LMFit fitter methods.

data_err: np.array(float)

An array containing the uncertainty on the data to be fitted. It is also stored as a one-dimensional array regardless of the type or dimensionality of the initial data.

noise: np.array(float) or None

If loaded, an array containing the background spectrum, including only the channels noticed in the fit.

noise_err: np.array(float or None)

If loaded, an array containing the sqrt of the background counts, only in the channels noticed during the fit. Used to compute the fit statistic.

_data_unmasked, _data_err_unmasked, _noise_unmasked: np.array(float)

The arrays of every data bin, its error and (if loaded) the backgruond, regardless of which ones are ignored or noticed during the fit. Used exclusively to enable book keeping internal to the fitter class.

Attributes inherited from EnergyDependentFit:

energs: np.array(float)

The array of physical photon energies over which the model is computed. Defined as the middle of each bin in the energy range stored in the instrument response provided.

energ_bounds: np.array(float)

The array of energy bin widths, for each bin over which the model is computed. Defined as the difference between the uppoer and lower bounds of the energy bins stored in the insrument response provided.

ear: np.array(float)

The array of energy bin bounds, for each bin over which the model is computed. Only necessary when calling Xspec models due to their unique input structure.

ebounds: np.array(float)

The array of energy channel bin centers for the instrument energy channels, as stored in the instrument response provided. Only contains the channels that are noticed during the fit.

ewidths: np.array(float)

The array of energy channel bin widths for the instrument energy channels, as stored in the instrument response provided. Only contains the channels that are noticed during the fit.

ebounds_mask: np.array(bool)

The array of instrument energy channels that are either ignored or noticed during the fit. A given channel i is noticed if ebounds_mask[i] is True, and ignored if it is false.

n_chans: int

The number of channels that are to be noticed during the fit.

_all_chans: int

The total number of channels in the loaded response matrix.

_emin_unmasked, _emax_unmasked, _ebounds_unmasked, _ewidths_unmasked: np.array(float)

The array of every lower bound, upper bound, channel center and channel widths stored in the response, regardless of which ones are ignored or noticed during the fit. Used exclusively to facilitate book-keeping internal to the fitter class.

Other attributes:

response: nDspec.ResponseMatrix

The instrument response matrix corresponding to the spectrum to be fitted. It is required to define the energy grids over which model and data are defined.

exposure: np.float

The exposure time of the observation. Only used for calculating Poisson-type likelihoods.

eval_model(params=None, ear=None, fold=True, mask=True)

This method is used to evaluate and return the model values for a given set of parameters, over a given model energy grid. By default it will evaluate the model over the energy grid defined in the response, using the parameters values stored internally in the model_params attribute, without folding the model through the response.

Parameters:

params: lmfit.Parameters, default None

The parameter values to use in evaluating the model. If none are provided, the model_params attribute is used.

ear: np.array(float), default None

The array of photon energy channel edges over which to evaluate the model. If none are provided, the same grid contained in the instrument response is used.

fold: bool, default True

A boolean switch to choose whether to fold the evaluated model through the instrument response or not. Not that in order for the model to be folded, the energy grid over which it is defined MUST be identical to that stored in the response matrix/class.

mask: bool, default True

A boolean switch to choose whether to mask the model output to only include the noticed energy channels, or to also return the ones that have been ignored by the users.

Returns:

model: np.array(float)

The model evaluated over the given energy grid, for the given input parameters.

plot_data(units='data', plot_bkg=False, return_plot=False)

This method plots the spectrum loaded by the user as a function of energy. It is possible to plot both in detector and ``unfolded’’ space, with the caveat that unfolding data is EXTREMELY dangerous and should be interpreted with care (or not at all).

The definition of unfolded data is subjective; nDspec adopts the same convention as ISIS, and defines an unfolded count spectrum Uf(h) as a function of energy channel h as : Uf(h) = C(h)/sum(R(E)), where C(h) is the detector space spectrum, R(E) is the instrument response and sum denotes the sum over energy bins. This definition has the advantage of being model-independent and is analogous to the Xspec (model-dependent) definition of unfolding data when the model is a constant.

It is also possible to return the figure object, for instance in order to save it to file.

Parameters:

units: str, default=”data”

The units to use for the y axis. units=”data”, the detector, plots the data in detector space in units of counts/s/keV. units=”unfold” instead plots unfolded data and follows the Xspec convention for the y axis - the y axis is in units of counts/s/keV/cm^2, times one additional factor “keV” for each “e” that appears in the string. For instance, units=”eeunfold” plots units of kev^2 counts/s/keV/cm^2, i.e. units of nuFnu.

plot_bkg; str, default=”False:

A boolean to choose whether you want to plot the background

return_plot: bool, default=False

A boolean to decide whether to return the figure objected containing the plot or not.

Returns:

fig: matplotlib.figure, optional

The plot object produced by the method.

plot_model(plot_data=True, plot_components=False, plot_bkg=False, params=None, units='data', residuals=None, return_plot=False)

This method plots the model defined by the user as a function of energy, as well as (optionally) its components, and the data plus model residuals. It is possible to plot both in detector and ``unfolded’’ space, with the caveat that unfolding data is EXTREMELY dangerous and should be interpreted with care (or not at all).

The definition of unfolded data is subjective; nDspec adopts the same convention as ISIS, and defines an unfolded count spectrum Uf(h) as a function of energy channel h as : Uf(h) = C(h)/sum(R(E)), where C(h) is the detector space spectrum, R(E) is the instrument response and sum denotes the sum over energy bins. This definition has the advantage of being model-independent and is analogous to the Xspec (model-dependent) definition of unfolding data when the model is a constant.

It is also possible to return the figure object, for instance in order to save it to file.

Parameters:

plot_data: bool, default=True

If true, both model and data are plotted; if false, just the model.

plot_components: bool, default=False

If true, the model components are overplotted; if false, they are not. Only additive model components will display their values correctly.

plot_bkg; str, default=False:

A boolean to choose whether you want to plot the background

params: lmfit.parameters, default=None

The parameters to be used to evaluate the model. If False, the set of parameters stored in the class is used

units: str, default=”data”

The units to use for the y axis. units=”data”, the detector, plots the data in detector space in units of counts/s/keV. units=”unfold” instead plots unfolded data and follows the Xspec convention for the y axis - the y axis is in units of counts/s/keV/cm^2, times one additional factor “keV” for each “e” that appears in the string. For instance, units=”eeunfold” plots units of kev^2 counts/s/keV/cm^2, i.e. units of nuFnu.

residuals: str, default=None

The units to use for the residuals. If residuals=”chisq”, the plot shows the residuals in units of data-model/error; if residuals=”ratio”, the plot instead uses units of data/model; if residuals “cstat”, the plot shows the contribution of each bin to the Cash statistic. If residual units are not specified, they are computed from the likelihood set by the user.

return_plot: bool, default=False

A boolean to decide whether to return the figure objected containing the plot or not.

Returns:

fig: matplotlib.figure, optional

The plot object produced by the method.

set_data(response, data, background=None)

This method sets the data to be fitted, its error, and the energy and channel grids given an input spectrum and its associated response matrix.

If the file provided was grouped with heatools, the method loads the grouped data and adjusts the channel grid automatically. The data is assumed to be background-subtracted (or to have negligible background).

Parameters:

response: nDspec.ResponseMatrix

An instrument response (including both rmf and arf) loaded into a nDspec ResponseMatrix object.

data: str

A string pointing to the path of an X-ray spectrum file, stored in a type 1 OGIP-formatted file (such as a pha file produced by a typical instrument reduction pipeline).

background: str, default None

a string pointing to the path of an X-ray spectrum background file, stored in a type 1 OGIP-formatted file. If not provided, the software assumes the data is either already background-subtracted, or that the user wants to ignore or model the background themselves.

set_fit_statistic(stat)

This method is used to set the statistic to be optimized during the fit. By default, the optimizer will optimize the chi-squared statistic.

Parameters:

stat: str

A string with the name of the fit statistic to be used. Supported statistics currently are “chisq” (the standard chi squared statistic, appropriate for data in the Gaussian regime) and “cstat” (the Cash statistic, see https://ui.adsabs.harvard.edu/abs/1979ApJ…228..939C/abstract, appropriate for Poisson-distributed data).

set_response(response)

This method sets the response matrix for the observation. It defines the energy grids over which model and data are defined. Generally, this method should only be called if the user is intending to simulate data from a model, as the response is not rebinned to reflect the data loaded by the user. Use the set_data method instead to set the data and response together.

Parameters:

response: nDspec.ResponseMatrix

An instrument response (including both rmf and arf) loaded into a nDspec ResponseMatrix object.

FitCrossSpectrum Class

class ndspec.FitCrossSpectrum.FitCrossSpectrum(likelihood='chisq')

Least-chi squares fitter class for the cross spectrum. The class supports both one-dimensional data between a reference and subject band as a function of Fourier frequency, and two-dimensional data between many subjects bands and a common reference band as a function of both energy and Fourier frequency.

Given an instrument response, an input spectrum, its error and a model, this class handles fitting internally using the lmfit library. The model can either be in the form of the actual values of the (complex) cross spectrum, a transfer function (defined in the Fourier domain), or an impulse response function (defined in the time domain). In all cases, the conversion from model units to cross spectral products (e.g. lag vs frequency, or real part vs energy in a fixed frequency interval) is handled automatically by the class.

Due to the variety of spectral-timing products related to or derived from the cross spectrum, users have several ways of defining the input data. In general, however, they are encouraged to do their own analysis separately using e.g. stingray, before moving on to fitting it.

model: lmfit.CompositeModel

A lmfit CompositeModel object, which contains a wrapper to the model component(s) one wants to fit to the data.

model_params: lmfit.Parameters

A lmfit Parameters object, which contains the parameters for the model components.

likelihood: str

A string that allows to switch between different fit statistics; which one is available depends on the type of fitter object. Uses chi-squared likelihood by default. Users can set different likelihoods either at initialization or with the appropriate setter method.

custom_likelihood: function

A function users can set to bypass the supported likelihoods and instead provide their own.

custom_args: tuple

A tuple including any custom arguments (in addition to the data and model values to be compared) necessary to calculate the custom likelihood

fit_result: lmfit.MinimizeResult

A lmfit MinimizeResult, which stores the result (including best-fitting parameter values, fit statistics etc) of a fit after it has been run.

data: np.array(float)

An array storing the data to be fitted. If the data is complex and/or multi-dimensional, it is flattened to a single dimension in order to be compatible with the LMFit fitter methods.

data_err: np.array(float)

An array containing the uncertainty on the data to be fitted. It is also stored as a one-dimensional array regardless of the type or dimensionality of the initial data.

_data_unmasked, _data_err_unmasked: np.array(float)

The arrays of every data bin and its error, regardless of which ones are ignored or noticed during the fit. Used exclusively to enable book keeping internal to the fitter class.

energs: np.array(float)

The array of physical photon energies over which the model is computed. Defined as the middle of each bin in the energy range stored in the instrument response provided.

energ_bounds: np.array(float)

The array of energy bin widths, for each bin over which the model is computed. Defined as the difference between the uppoer and lower bounds of the energy bins stored in the insrument response provided.

ear: np.array(float)

The array of energy bin bounds, for each bin over which the model is computed. Only necessary when calling Xspec models due to their unique input structure.

ebounds: np.array(float)

The array of energy channel bin centers for the instrument energy channels, as stored in the instrument response provided. Only contains the channels that are noticed during the fit.

ewidths: np.array(float)

The array of energy channel bin widths for the instrument energy channels, as stored in the instrument response provided. Only contains the channels that are noticed during the fit.

ebounds_mask: np.array(bool)

The array of instrument energy channels that are either ignored or noticed during the fit. A given channel i is noticed if ebounds_mask[i] is True, and ignored if it is false.

n_chans: int

The number of channels that are to be noticed during the fit.

_all_chans: int

The total number of channels in the loaded response matrix.

n_bins: int

The number of noticed channels times the number of noticed bins in Fourier frequency.

_all_bins: int

The total number of channels, times the total number of bins in Fourier frequency.

_emin_unmasked, _emax_unmasked, _ebounds_unmasked, _ewidths_unmasked: np.array(float)

The array of every lower bound, upper bound, channel center and channel widths stored in the response, regardless of which ones are ignored or noticed during the fit. Used exclusively to facilitate book-keeping internal to the fitter class.

_freqs_unmasked: np.array(float)

If the data and model explicitely depend on Fourier frequency (e.g. a power spectrum), this is the array of Fourier frequency over which all data and model are defined, including bins that are ignored in the fit.

If instead the data depends from some other energy (e.g. energy), it contains both noticed and ignored frequency intervals over which to produce spectral-timing products. For example, a user might input a set of 7 ranges of frequencies to calculate lag energy spectra, but only want to consider the first and last 3, and ignore the middle one.

freqs_mask np.array(bool)

The array of Fourier frequencies that are either ignored or noticed during the fit. A given channel i is noticed if freqs_mask[i] is True, and ignored if it is false.

n_freqs: int

The number of Fourier frequency bins that are noticed in the fit.

response: nDspec.ResponseMatrix

The instrument response matrix corresponding to the spectrum to be fitted. It is required to define the energy grids over which model and data are defined.

units: str

A string that checks the units which the user is providing - “lags” for fitting lag spectra alone, “polar” or fitting modulus and phase together, and “cartesian” for fitting real and imaginary parts together.

ref_band: [np.float,np.float]

The minimum/maximum energy bounds over which to take the reference band. Necessary to calculate spectral timing products (like lag spectra) from the input model.

freqs: np.array(float)

If the data explicitely depends on Fourier frequency, it is the range of Fourier frequencies over which both the data and model are defined. Otherwise, it is the internal Fourier frequency grid over which the model is computed before being converted into spectral-timing products (e.g. lag-energy spectra).

freq_bounds: np.array(float)

The array of Fourier frequency bounds over which energy-dependent data and model are averaged over, in order to handle energy-dependent spectral-timing products.

_times: np.array(float)

The array of times corresponding to the Fourier frequency array. Used internally for model calculations/book-keeping.

crossspec: nDspec.CrossSpectrum

A nDspec CrossSpectrum object used to store model evaluations, and to convert model evaluations into spectral-timing data products.

renorm_phase: bool

Allows users to apply a small phase renormalization when fitting energy dependent products. This is necessary to account for imperfections in the instrument response/calibration. For more discussion, see Appendix E in Mastroserio et al. 2018: https://ui.adsabs.harvard.edu/abs/2018MNRAS.475.4027M/abstract This setting will NOT affect the modulus of a cross spectrum, only the phase (and therefore it will affect the real and imaginary parts).

renorm_modulus: bool

Allows users to apply a small modulus renormalization when fitting energy dependent products. This can be useful when defining models from transfer or impulse response functions, in order to re-normalize the model in each Fourier frequency bin to match the data. Physically, this allows one to take into account differences between the power spectrum shape assumed in the model to calculate the spectral timing products (e.g. modulus vs energy), and the ``true’’ underyling power spectrum in the source. This setting will NOT affect the phase of a cross spectrum, only the modulus (and therefore it will affect the real and imaginary parts).

dependence: str

A string recording whether the data and model to be fitted are in units of frequency or energy (e.g. lag frequency vs lag energy spectra).

_supported_coordinates: str

A string that checks the units models/data can be defined as. “lags” is for fitting lag spectra alone, “polar” is for fitting modulus and phase together, and “cartesian” is for fitting real and imaginary parts together.

_supported_models: str

A string that checks the type of model defined. “cross” indicates models that already return the full cross spectrum, and thus need limited operations applied before being compared to the data. “transfer” indicates models of transfer functions, which need to be crossed with a reference band to calculate the cross spectrum. “irf” indicates a models of impulse response functions, which need to be Fourier transformed and then crossed with a reference band.

_supported_products: str

A string that checks whether the data provided (e.g. lags) is a function of Fourier frequency or energy.

needbkg: bool, default=True

A boolean that indicates whether the background is needed for the simulation of the cross spectrum. If it is set to True, the class will attempt to read the background from a file specified by the user and set this attribute to False. If it is set to False, it is assumed that the background file has already been read in and will not be read again.

eval_model(params=None, fold=True, mask=True)

This method is used to evaluate and return the model values for a given set of parameters, over the internal energy and frequency grids. By default the model is evaluated using the parameters values stored internally in the model_params attribute. The model is always folded through the instrument response, returning either all or only the noticed channels. The reference band is always that set in set_data.

Parameters:

params: lmfit.Parameters, default None

The parameter values to use in evaluating the model. If none are provided, the model_params attribute is used.

fold: bool, default True

A boolean switch to choose whether to fold the model through the instrument response or not.

mask: bool, default True

A boolean switch to choose whether to mask the model output to only include the noticed energy channels, or to also return the ones that have been ignored by the users.

Returns:

model: np.array(float)

The model evaluated over the given energy grid, for the given input parameters.

plot_data_1d(return_plot=False)

This method plots the cross spectrum loaded by the user as a function of the unit dependence specified (ie, Fourier frequency or energy). If the data is loaded is two-dimensional - as is the case for frequency dependent products in multiple subject bands, or energy dependent products in multiple Fourier frequency bins, the method plots every loaded spectrum in a single plot.

Regardless of the data format, only the noticed energy channels are plotted. Note that if the user is ignoring energy bins that are not at the limit of the channel grid (e.g., between 5 and 7 keV for a grid of channels between 0.5 and 10 keV), then the automated legend will not label the spectra correctly.

It is also possible to return the figure object, for instance in order to save it to file.

Parameters:

return_plot: bool, default=False

A boolean to decide whether to return the figure objected containing the plot or not.

Returns:

fig: matplotlib.figure, optional

The plot object produced by the method.

plot_data_2d(use_phase=False, return_plot=False)

This method plots the cross spectrum loaded by the user in two dimensions as a function of both Fourier frequency or energy. Regardless of the data format, only the noticed energy channels are plotted. Note that due to how matplotlib handles two-dimensional plots, the bounds of the bins shown in the plot will differ slightly from those where the data is defined.

It is also possible to return the figure object, for instance in order to save it to file.

Parameters:

use_phase: bool, default=False

A boolean used exclusively when plotting time lags. If it is true, it converts the lags to phases for ease of visualization over a large range in lag timescales.

return_plot: bool, default=False

A boolean to decide whether to return the figure objected containing the plot or not.

Returns:

fig: matplotlib.figure, optional

The plot object produced by the method.

plot_model_1d(plot_data=True, params=None, residuals='chisq', return_plot=False)

This method plots the model defined by the user as a function of the unit dependence specified (ie, Fourier frequency or energy).

By default the method also plots the data loaded as well as the residuals. Additionally, by default the model is evaluated using the parameter values currently stored in the class, but it is also posssible to pass a different set of parameters for model evaluation.

If the data is loaded is two-dimensional - as is the case for frequency dependent products in multiple subject bands, or energy dependent products in multiple Fourier frequency bins, the method plots every loaded spectrum in a single plot.

Regardless of the data format, only the noticed energy channels are plotted. Note that if the user is ignoring energy bins that are not at the limit of the channel grid (e.g., between 5 and 7 keV for a grid of channels between 0.5 and 10 keV), then the automated legend will not label the spectra correctly.

It is also possible to return the figure object, for instance in order to save it to file.

Parameters:

plot_data: bool, default=True

If true, both model and data are plotted; if false, just the model.

params: lmfit.parameters, default=None

The parameters to be used to evaluate the model. If False, the set of parameters stored in the class is used

residuals: str, default=”chisq”

The units to use for the residuals. If residuals=”chisq”, the plot shows the residuals in units of data-model/error; if residuals=”ratio”, the plot instead uses units of data/model. If using a custom likelihood, the residuals are computed from it.

return_plot: bool, default=False

A boolean to decide whether to return the figure objected containing the plot or not.

Returns:

fig: matplotlib.figure, optional

The plot object produced by the method.

plot_model_2d(params=None, use_phase=False, residuals='chisq', return_plot=False)

This method plots the model and data loaded by the user in two dimensions as a function of both Fourier frequency or energy. Regardless of the data format, only the noticed energy channels are plotted. Note that due to how matplotlib handles two-dimensional plots, the bounds of the bins shown in the plot will differ slightly from those where the data is defined.

It is also possible to return the figure object, for instance in order to save it to file.

Parameters:

params: lmfit.parameters, default=None

The parameters to be used to evaluate the model. If False, the set of parameters stored in the class is used

use_phase: bool, default=False

A boolean used exclusively when plotting time lags. If it is true, it converts the lags to phases for ease of visualization over a large range in lag timescales.

residuals: str, default=”chisq”

The units to use for the residuals. If residuals=”chisq”, the plot shows the residuals in units of data-model/error; if residuals=”ratio”, the plot instead uses units of data/model. If using a custom likelihood, the residuals are computed from it.

return_plot: bool, default=False

A boolean to decide whether to return the figure objected containing the plot or not.

Returns:

fig: matplotlib.figure, optional

The plot object produced by the method.

renorm_mods(switch)

Setter method to enable the modulus renormalization when fitting energy depenent products. This renormalization is intended to correct for not knowing the correct shape of the power spectrum responsible the observed variability. See Mastroserio et al. 2018 for further discussion: https://ui.adsabs.harvard.edu/abs/2018MNRAS.475.4027M/abstract This setting will NOT affect the phase of a cross spectrum, only the modulus (and therefore it will affect the real and imaginary parts).

Parameters:

switch: bool

A boolean to track whether modulus renormalization is enabled or not. If it is, the method modifies the defined model and its parameters automatically.

renorm_phases(switch)

Setter method to enable the phase renormalization when fitting energy depenent products. This renormalization is intended to correct for small uncertainties in the instrument response function, which can affect the phase of the cross spectrum. For more discussion, see Appendix E in Mastroserio et al. 2018: https://ui.adsabs.harvard.edu/abs/2018MNRAS.475.4027M/abstract This setting will NOT affect the modulus of a cross spectrum, only the phase (and therefore it will affect the real and imaginary parts).

Parameters:

switch: bool

A boolean to track whether phase renormalization is enabled or not. If it is, the method modifies the defined model and its parameters automatically.

set_background(bkg_file_path)

This method is used to set the background file to be used in the simulation of lag spectra. The background file should contain the background counts in each energy channel.

Parameters:

bkg_file: str

The path to the background file containing the background counts for each energy channel.

set_coordinates(units)

This method allows users to specify the coordinate system of the data they intend to fit. The possible choices are “lags”, for time lags alone, “polar”, for modulus and phase together, or “cartesian”, for real and imaginary parts together.

Parameters:

units: str

A string containing the dependence of the data. Must be one of the supported units stored in the class - ie “lags”, “polar” or “cartesian”.

set_data(response, ref_bounds, sub_bounds, data, data_err=None, freq_grid=None, time_grid=None, freq_bins=None, time_res=None, seg_size=None, norm=None)

This method is used to set the cross-spectrum data to be fitted. The exact data is determined by the set_product_dependence and set_coordinates setter methods.

The method requires an input instrument response, bounds defined for the reference band as well as the subject band(s), and the actual data. This can be in the form of an array, in which case users also need to specify the errors, or (only for frequency-dependent data) a stingray.events object. Additionally, users need to specify the time resolution and segment size of the lightcurves used to build the data.

When loading energy-dependent products (e.g. many lag-energy spectra), it is necessary to specify the Fourier frequency intervals over which each lag spectrum is computed, using the freq_bins argument.

When loading freqency-dependent products from a stingray event file, it is possible to specify the normalization to be used. By default, this will be asbsolute rms normalization.

Finally, for both products, users can specify by hand the Fourier frequency and time grids to be used internally for model computations. In some cases, this can be useful in speeding up model evaluations.

Parameters:

response: nDspec.ResponseMatrix

The instrument response matrix corresponding to the spectrum to be fitted. It is required to define the energy grids over which model and data are defined. It is rebinned automatically such that the subject/reference bands are in separate channels.

ref_bounds: [np.float,np.float]

The minimum/maximum energy bounds over which to take the reference band.

sub_bounds: np.array([float,float])

An array of minimum/maximum energy bounds over which each channel of interest is taken.

data: np.array(float) or stingray.events

The data to be fitted, either in the form of a one-dimensional array or a stingray event file. If passing an array, the data has to be stored such that all the real values or moduli are contained in the first half of the array, with the imaginary values or phases in the second half.

data_err: np.array(float), optional

Only required when not passing a stingray event file. Needs to be in the same format as the data array - all moduli/real parts first, all phases/imaginary parts second.

freq_grid: np.array(float), optional

The grid of Fourier frequency over which to compute the model. If it is not passed explicitely, it is computed from the time resolution and segment size arguments.

time_grid: np.array(float), optional

The grid of times used to generate the Fourier frequency grid. Used for internal model calcluations, if using an impulse response function and the users wishes to use the sinc decomposition method.

freq_bins: np.array(float), optional

The Fourier frequency grids over which energy-dependent data has been averaged. Required for fitting energy-dependent data (e.g. lag energy spectra), and not used for frequency-dependent data.

time_res: np.float, optional

The time resolution of the lightcurves used to build the data. It is necessary to build the grid of Fourier frequency over which to compute the model.

seg_size: np.float, optional

The size of the segments in which the lightcurves were divided to build the data. It is necessary to build the grid of Fourier frequency over which to compute the model.

norm; str, optionalm default = “abs”

The normalization of the data products, if they are calculated from a stingray event file. If not specified, absolute rms normalization is used.

set_model(model, model_type='irf', params=None)

This method is used to pass the model users want to fit to the data. Users also need to specify what quantity the model actually computes - a time-dependent impulse respnse function, a Fourier-frequency dependent function, or an explicit value for the cross spectrum. Based on the model type, the class then converts model output to the appropriate spectral timing products automatically. Optionally it is also possible to pass the initial parameter values of the model.

Parameters:

model: lmfit.CompositeModel

The lmfit wrapper of the model one wants to fit to the data.

model_type: str, default=”irf”

A string describing the type of models users want to define. Options are “irf” for an impuplse response function, “transfer” for a transfer functino, and “cross” for the actual cross spectrum.

params: lmfit.Parameters, default: None

The parameter values from which to start evalauting the model during the fit. If it is not provided, all model parameters will default to 0, set to be free, and have no minimum or maximum bound.

set_product_dependence(depend)

This method allows users to specify whether the data they intend to fit is a function of Fourier frequency (e.g. lag frequency spectra), or energy (e.g. lag vs energy spectra). Polarimetric data is not supported at this time.

Parameters:

depend: str

A string containing the dependence of the data. Must be one of the supported dependence stored in the class - ie “frequency” or “energy”.

set_psd_weights(psd_weights)

This method is necessary when users define models from an impulse response or transfer functins, and sets the power spectrum used as weights when calculating spectral timing products.

Parameters:

input_power: np.array(float) or PowerSpectrum

Either an array of size (len(freqs)) that is to be used as the weighing power spectrum when computing the cross spectrum, or an nDspec PowerSpectrum object. Both have to be defined over the class internal Fourier frequency grid.

FitTwoD Class

class ndspec.FitTwoD.FitTwoD(likelihood='chisq')

This class is designed for fitting generic types of two-dimensional data, regardless of what units it may be in. Models used in this fitter are expected to be provided already in the same unit as the data. Common examples of using this class might be time-dependent spectroscopy, or fitting a dynamical power spectrum.

As an exception, users can optionally pass an instrument response matrix object, in which case the y axis is assumed to be in units of photon channels and the model should produce units of integrated photon flux over that axis.

Attributes inherited from SimpleFit:

model: lmfit.CompositeModel

A lmfit CompositeModel object, which contains a wrapper to the model component(s) one wants to fit to the data.

model_params: lmfit.Parameters

A lmfit Parameters object, which contains the parameters for the model components.

likelihood: str

A string that allows to switch between different fit statistics; which one is available depends on the type of fitter object. Uses chi-squared likelihood by default. Users can set different likelihoods either at initialization or with the appropriate setter method.

custom_likelihood: function

A function users can set to bypass the supported likelihoods and instead provide their own.

custom_args: tuple

A tuple including any custom arguments (in addition to the data and model values to be compared) necessary to calculate the custom likelihood

fit_result: lmfit.MinimizeResult

A lmfit MinimizeResult, which stores the result (including best-fitting parameter values, fit statistics etc) of a fit after it has been run.

data: np.array(float)

An array storing the data to be fitted. Only contains noticed bins.

data_err: np.array(float)

An array containing the uncertainty on the data to be fitted. Only contains noticed bins.

_data_unmasked, _data_err_unmasked: np.array(float)

The arrays of every data bin and its error, regardless of which ones are ignored or noticed during the fit. Used exclusively to enable book keeping internal to the fitter class.

Other attributes:

rows, columns: int

Two integers keeping track of the number of rows and columns in the data loaded. Only tracks the number of rows/columns that are noticed during the fit.

column_grid: np.array(float)

The grid of values along which the x-axis of the model is computed.

column_mask: np.array(bool)

A masking array used to ignore or notice data bins along the x axis.

row_grid: np.array(float)

The grid of values along which the y-axis of the model is computed. If an instrument response is loaded, it contains the photon channels.

column_mask: np.array(bool)

A masking array used to ignore or notice data bins along the y axis.

_column_grid_unmasked,_row_grid_unmasked: np.array(float)

The unmasked arrays over which the entire data is define, whether they are noticed duing the fit or not.

_all_rows,_all_columns,_all_bins: int

Integers keeping track of all the rows, columns, or total number of data points, whether they are noticed in the fit or not.

dependence, units: str

Two strings used to specify the units the data is in, used to handle filtering in/out data bins in the fitter.

response: nDspec.ResponseMatrix, default None

The instrument response matrix corresponding to the data to be fitted. It is required to define the energy grids over which model and data are defined. The following arrays are only defined when a response matrix is loaded:

energs: np.array(float)

The array of physical photon energies over which the model is computed. Defined as the middle of each bin in the energy range stored in the instrument response provided.

energ_bounds: np.array(float)

The array of energy bin widths, for each bin over which the model is computed. Defined as the difference between the uppoer and lower bounds of the energy bins stored in the insrument response provided.

ear: np.array(float)

The array of energy bin bounds, for each bin over which the model is computed. Only necessary when calling Xspec models due to their unique input structure.

ebounds: np.array(float)

The array of energy channel bin centers for the instrument energy channels, as stored in the instrument response provided. Only contains the channels that are noticed during the fit.

ewidths: np.array(float)

The array of energy channel bin widths for the instrument energy channels, as stored in the instrument response provided. Only contains the channels that are noticed during the fit.

eval_model(params=None, column_grid=None, row_grid=None, fold=True, mask=True)

This method is used to evaluate and return the model values for a given set of parameters, over given row and column grids. If a response is loaded, by default it will evaluate the model over the energy grid defined in the response, using the parameters values stored internally in the model_params attribute, without folding the model through the response.

Parameters:

params: lmfit.Parameters, default None

The parameter values to use in evaluating the model. If none are provided, the model_params attribute is used.

column_grid: np.array(float), default None

The array of values in the x-axis over which to evaluate the model. If none are provided, the same grid contained in the fitter object is used.

row_grid: np.array(float), default None

The array of values in the y-axis over which to evaluate the model. If none are provided, the same grid contained in the fitter object is used. If the fitter contains an instrument response, this array is the energy grid contained in that response.

fold: bool, default True

A boolean switch to choose whether to fold the evaluated model through the instrument response or not. Not that in order for the model to be folded, the energy grid over which it is defined MUST be identical to that stored in the response matrix/class.

mask: bool, default True

A boolean switch to choose whether to mask the model output to only include the noticed energy channels, or to also return the ones that have been ignored by the users.

Returns:

model: np.array(float)

The model evaluated over the given energy grid, for the given input parameters.

ignore_columns(bound_lo, bound_hi)

This method adjusts the arrays stored such that they (and the fit) ignore selected columns based on their bounds.

Parameters:

bound_lofloat

Lower bound of ignored column interval.

bound_hifloat

Higher bound of ignored column interval.

ignore_rows(bound_lo, bound_hi)

This method adjusts the arrays stored such that they (and the fit) ignore selected columns rows on their bounds.

Parameters:

bound_lofloat

Lower bound of ignored rows interval.

bound_hifloat

Higher bound of ignored rows interval.

notice_columns(bound_lo, bound_hi)

This method adjusts the arrays stored such that they (and the fit) notice selected columns based on their bounds.

Parameters:

bound_lofloat

Lower bound of noticed column interval.

bound_hifloat

Higher bound of noticed column interval.

notice_rows(bound_lo, bound_hi)

This method adjusts the arrays stored such that they (and the fit) notice selected columns rows on their bounds.

Parameters:

bound_lofloat

Lower bound of noticed rows interval.

bound_hifloat

Higher bound of noticed rows interval.

plot_data(x_name=None, y_name=None, return_plot=False)

This method creates a 2d plot of the data loaded by the user as a function of whatever units are in used for the rows and columns.

It is also possible to return the figure object, for instance in order to save it to file.

Parameters:

x_name: str, default=”None”

The units to use to label the X axis.

y_name: str, default=”None”

The units to use to label the Y axis.

return_plot: bool, default=False

A boolean to decide whether to return the figure objected containing the plot or not.

Returns:

fig: matplotlib.figure, optional

The plot object produced by the method.

plot_model(residuals='chisq', params=None, x_name=None, y_name=None, return_plot=False)

This method creates a 2d plot of the data loaded by the user, the model for the given parameters (either passed as an object, or already loaded in the fitter), and residuals, as a unction of whatever units are in used for the rows and columns.

It is also possible to return the figure object, for instance in order to save it to file.

Parameters:

residuals: str, default=”chisq”

The units to be used in the residuals.

params: lmfit.parameters, default=None

The parameters to be used to evaluate the model. If False, the set of parameters stored in the class is used.

x_name: str, default=”None”

The units to use to label the X axis.

y_name: str, default=”None”

The units to use to label the Y axis.

return_plot: bool, default=False

A boolean to decide whether to return the figure objected containing the plot or not.

Returns:

fig: matplotlib.figure, optional

The plot object produced by the method.

set_data(data, data_err, column_grid, row_grid, response=None, noise=None, noise_err=None)

This method is used to pass the data users want to fit. The input has to be in the form of numpy arrays for the data, its errors, and the grids over which it is defined. Users can optionally pass an instrument response, as well as background arrays.

Parameters:

data, data_err: np.array([float,float])

Two-dimensional arrays containing the data to be fitted and its uncertainties.

column_grid, row_grid: np.array(float)

The arrays over which the data is defined. If users pass an instrument response (see below), then row_grid is the array of energy bin bounds, for each channel over which the data is defined. This is identical to the ‘ear’ array in Xspec models.

response: nDspec.ResponseMatrix, default None

An instrument response (including both rmf and arf) loaded into a nDspec ResponseMatrix object.

noise, noise_err: np.array(float), default None

Optional arrays including the background/noise floor and its error.

set_fit_statistic(stat)

This method is used to set the statistic to be optimized during the fit. By default, the optimizer will optimize the chi-squared statistic.

Parameters:

stat: str

A string with the name of the fit statistic to be used. Supported statistics currently are “chisq” (the standard chi squared statistic, appropriate for data in the Gaussian regime) and “cstat” (the Cash statistic, see https://ui.adsabs.harvard.edu/abs/1979ApJ…228..939C/abstract, appropriate for Poisson-distributed data).

JointFit Class

class ndspec.JointFit.JointFit

Generic joint fitting class. Use this class if you have multiple datasets that you want to fit jointly in some way.

Users can add other Fit…Objs from ndspec to an instance of this class, and that instance this will handle evaluating the model, sharing parameters and/or whole models between Fit….Objs objects, and perform inference and/or optimization. There is no restriction on the type or number of Fit…Objs that can be added.

Attributes:

jointdict{Fit… objects and/or list(Fit… objects)}

Dictionary containing named Fit… objects to be joint-fitted.

joint_params: dict{lists(str)}

Dictionary containing the names of model parameters for each distinct fit object.

fit_result: lmfit.MinimizeResult

A lmfit MinimizeResult, which stores the result (including best-fitting parameter values, fit statistics etc) of a fit after it has been run.

model_params: lmfit.Parameters

The parameters of all the fitter objects combined, which are used during the joint fit.

energy_grid: dict{np.arrays}

An optional dictionary containing the details of a user-defined grid of energy bins and bounds over which to evaluate a model. Usable only with time-averaged spectra using a single shared model. Using a shared grid can improve performance in some cases.

shared_energy_grid: bool

A boolean that tracks whether the joint fitter object has a shared energy grid for evaluating the time-averaged spectrum model or not.

shared_model: lmfit.CompositeModel

The model shared among all FitTimeAvgSpectrum objects to be evaluated on the same energy grid.

spec_renorm_model: lmfit.Model

An optional lmfit model object used to include a constant component in a set of chosen fitters for time-averaged spectra. Used to account for instrument flux cross-calibration.

shared_keys: list

A list containing the keys (names) of the parameters that are identical and therefore shared between the loaded fitter objects.

renorm_names: list[str]

A list with the names of the time-averaged fitter objects to which the cross-calibration constants should be applied.

add_fitobj(fitobj, name)

Adds one or more fitters to the joint fit.

Parameters

fitobjFit… object or list of Fit… objects

the Fit… object(s) being included in the joint fit.

name: str

name(s) to be assigned to the fitter objects loaded. These will be used as keys of a dictionary to retrieve the individual fitter objects and their methods.

all_plots(units, residuals='chisq', plot_bkg=None, return_plot=False)

This method loops over all stored fitter objects and plots the data, model (given the parameters stored), and residuals for all the fits separately. For two-dimensional fits (like a cross spectrum), the method still plots one-dimensional plots. This method is meant for quick-look analysis of a joint fit.

Parameters:

units: str

The units to use for the y axis. For more info, see the documentation of the individual fitter classes. Has no impact if the fitter being looked over is a cross spectrum.

residuals: str, default=”chisq”

The units to use for the residuals. If residuals=”chisq”, the plot shows the residuals in units of data-model/error; if residuals=”ratio”, the plot instead uses units of data/model. For cross spectra this key word is ignored and only delta chi residuals can be shown.

plot_bkg; str, default=False:

A boolean to choose whether you want to plot the background

return_plot: bool, default=False

A boolean to decide whether to return the figure objected containing the plot or not.

eval_model(params=None, names=None, flatten=True)

This method is used to evaluate and return the model values of models stored in the fitter objects that make up the joint fit.

Parameters:

params: lmfit.Parameters, default None

The parameters to use in evaluating the model/models. If none are provided, the model_params attribute is used.

names: list(str), default None

Names of the fitters that should evalualate their models. Defaults to evaluating the models of all fitters stored in the joint fit.

flatten: bool, default True

A boolean to switch between returning model evaluations as a dictionary or numpy array (see below).

Returns:

model_hierarchy: either dict(np.array(float)) or np.array(float)

Models are evaluated and returned either as a dictionary, with keys defined by the fitter names, or by flattened numpy arrays. The former allows easy access to the evaluated model for each fitter, the latter is necessary for lmfit optimzers and/or likelihood calculations.

fit_data(algorithm='leastsq', names=None, report_result=True)

This method attempts to minimize the residuals of the model with respect to the data defined by the user. The fit either starts from the set of parameters defined with .set_params(), or from the parameters set in the individual fitter obejcts loaded. Once the algorithm has completed its run, it optionally prints to terminal the best-fitting parameters, fit statistics, and parameter values.

Parameters:

algorithm: str, default=”leastsq”

The fitting algorithm to be used in the minimization. The possible choices are detailed on the LMFit documentation page: https://lmfit.github.io/lmfit-py/fitting.html#fit-methods-table.

names: list(str), default None

names of the fit objects that should be part of the fit. Defaults to using all the fitters stored in the JointFit object instance.

get_residuals(params=None, names=None, flatten=True)

This methods return the residuals based on the statistic, model and data defined in each of the loaded fitter objects.

Parameters:

params: lmfit.Parameters or list[lmfit.Parameters]

The parameter values to use in evaluating the models. Defaults to those stored in the JointFit instance.

names: list(str), default None

names of the models that should be evalualated. Defaults to evaluating all models.

flatten: bool, default True

A boolean to switch between returning model evaluations as a dictionary or numpy array (see below).

Returns:

residuals: np.array(float) or dict{np.array(float)}

Either an array of the same size as the total number of data points, or a dictionary with keys named after each loaded SimpleFit instance, containing the model residuals in each channel.

joint_plot(units, residuals='chisq', plot_bkg=False, xrange=None, yrange=None, return_plot=False, names=None)

This method loops over all stored fitter objects and plots the data, model (given the parameters stored), and residuals for all the fits together. Note that this is useful only if the data you are trying to plot is of the same time (e.g. all time-averaged spectra).

Parameters:

units: str

The units to use for the y axis. For more info, see the documentation of the individual fitter classes.

residuals: str, default=”chisq”

The units to use for the residuals. If residuals=”chisq”, the plot shows the residuals in units of data-model/error; if residuals=”ratio”, the plot instead uses units of data/model. For cross spectra this key word is ignored and only delta chi residuals can be shown.

plot_bkg; str, default=False:

A boolean to choose whether you want to plot the background

xrange, yrange: (float, float)

The limits of the plot on the x and y axis

return_plot: bool, default=False

A boolean to decide whether to return the figure objected containing the plot or not.

names: list(str)

A list of the fitters to be included in the joint plot. By default, all the loaded fitters are included.

print_fit_report()

This method prints the current fit result.

print_models(names=None)

Prints the models contained within the joint fit. Defaults to printing all models, but users can filter models using the names parameter.

Parameters

namesstr or list(str), optional

names of the models to be printed. The default is to print all models.

renorm_timeavg(switch, names=None)

Setter method to enable the fitter objects passed with the “names” attribute to include an additional constant component, in other to account e.g. for the cross-calibration uncertainty between instruments. Only applicable to time-averaged spectra.

Parameters:

switch: bool

A boolean to track whether the spectra renormalization is enabled or not. If it is, the method modifies the defined model and its parameters automatically.

name: list(str) or str, default None

A list of strings with the name of the fitter objects to which users wish to apply a cross calibration constant. By default this is None and all the spectra receive the constant.

set_energy_grid(grid_bounds)

This method defines a custom energy grid over which to evaluate all loaded time-averaged spectra. The new grid MUST cover a wider energy range than all the ones in the time-averaged spectra loaded in the joint fitter instance.

Parameters:

grid_bounds: np.array(float)

An array of energy bounds, starting from the bottom edge of the first bin and finishing up to the top edge of the last bin in the new grid.

set_params(params)

This method is used to set the model parameter names and values. It can be used both to initialize a fit, and to test different parameter values before actually running the minimization algorithm.

Parameters:

params: lmfit.parameter

The parameter values from which to start evalauting the model during the fit.

Sampling functions

ndspec.SamplingUtils.set_sampling_priors(fitobj, priors)

This function is used to set the priors to be used with emcee sampling. These priors are saved in a global variable called sampling_priors; therefore, users should never re-use the variable name sampling_priors in their code.

Parameters:

fitobj: ndspec.Fit…Object or ndspec.JointFit

Object containing the data, specified model, and parameters

priors: dict

A dictionary of priors to be used in emcee. The key of each dictionary should be the name of the parameter. Each key should contain an object with a method called “logprob”, which returns the (negative) logarithm of the prior evaluated at a given point.

ndspec.SamplingUtils.set_sampling_model(fitobj)

This function is used to set the model to be used with emcee sampling. This model is saved in a global variable called sampling_model; therefore, users should never re-use the variable name sampling_model in their code.

Parameters:

fitobj: ndspec.Fit…Object or ndspec.JointFit

Object containing the data, specified model, and parameters

ndspec.SamplingUtils.set_sampling_data(fitobj)

This function is used to set the data and its error to be used with emcee sampling. These are saved in global variables called sampling_data and sampling_data_err; therefore, users should never re-use the variable names sampling_data and sampling_data_err in their code. If the fitter object includes noise (e.g. a background spectrum) and exposure times, these are included as well.

Parameters:

fitobj: ndspec.Fit…Object or ndspec.JointFit

Object containing the data, specified model, and parameters

ndspec.SamplingUtils.set_sampling_parameters(params)

This function is used to set the parameters of the model to be used with emcee sampling. The parameter object (containing all parameters), as well as the names and values of the variable parameters are saved in global variables called sampling_names, sampling_values and sampling_params; therefore, users should never re-use these variable names in their code.

Parameters:

params: lmfit.Parameters

The lmfit parameters object used in the model, including those kept constant.

Returns:

theta: np.array

A numpy array containing the values of the free parameters in the model.

ndspec.SamplingUtils.initialise_mcmc(fitobj, priors)

This function is used to initialise an MCMC run. The Fit…Object can be any of the particular data products, or a JointFit object containing multiple FitObjects.

Parameters:

fitobj: ndspec.Fit…Object or ndspec.JointFit

Object containing the data, specified model, and parameters

priors: dict

A dictionary of priors to be used in emcee. The key of each dictionary should be the name of the parameter. Each key should contain an object with a method called “logprob”, which returns the (negative) logarithm of the prior evaluated at a given point.

Returns:

theta: np.array

A numpy array containing the values of the free parameters in the model.

ndspec.SamplingUtils.reflect_parameter(x, a, b)

This function is used to bounce a parameter value off a hard-limit, defined by e..g the limits of a uniform prior. It is useful to improve the behavior and acceptance rate of a mcmc chain

Parameters:

x: np.array(float) or (float)

A value, or array of values, for a parameter which may or may not exceed the limits defined by the priors

a: np.array(float) or (float)

A value, or array of values, containing the lower limit allowed for each parameter

b: np.array(float) or (float)

A value, or array of values, containing the upper limit allowed for each parameter

Returns:

y: np.array(float) or (float)

The parameter value to be used in the model, “bounced” of the hard limit if it exceeds it - e.g. if x=1.2, a=0, b=1, then y = 0.8.

class ndspec.SamplingUtils.priorUniform(min, max)

This class is used to compute a uniform prior distribution during Bayesian sampling, for a given model parameter.

Parameters:

min: float

The lower bound of the distribution.

max: float

The upper bound of the distribution.

logprob(theta)

This method returns the log probability of the distribution - in this case, an (arbitrary, for the purpose of likelihood optimization) constant.

Parameters:

theta: float

The parameter value for which the likelihood is to be computed.

Returns:

The value of the likelihood for the input parameter.

transform_prior(quantile)

This method returns the parameter value defined from the quantile of the cumulative distribution defined in the prior. It effectively converts an interval from 0-1 to one that samples the prior distribution and is used by nested sampling to define prior likelihoods.

Parameters:

quantile: float, 0-1

The quantile for which the prior is to be computed

Returns:

prior: float

The parameter value for the chosen parameter prior distribution.

class ndspec.SamplingUtils.priorLogUniform(min, max)

This class is used to compute a log-uniform prior distribution (in base e) during Bayesian sampling, for a given model parameter.

Parameters:

min: float

The lower bound of the distribution

max: float

The upper bound of the distribution

logprob(theta)

This method returns the log probability of the distribution - in this case, an (arbitrary, for the purpose of likelihood optimization) constant. More explicitely, if x is our parameter and log10(x) is uniform, then p(log10(x)) = const, p(x) = p(log10(x))*dlog10(x)/dx = const/x. Therefore, the log-probability is (minus a constant) log10(p(x)) = log10(1/x) = -log10(x).

Parameters:

theta: float

The parameter value for which the likelihood is to be computed.

Returns:

The value of the likelihood for the input parameter.

transform_prior(quantile)

This method returns the parameter value defined from the quantile of the cumulative distribution defined in the prior. It effectively converts an interval from 0-1 to one that samples the prior distribution and is used by nested sampling to define prior likelihoods.

Parameters:

quantile: float, 0-1

The quantile for which the prior is to be computed

Returns:

prior: float

The parameter value for the chosen parameter prior distribution.

class ndspec.SamplingUtils.priorNormal(mu, sigma)

This class is used to compute a normal prior distribution during Bayesian sampling, for a given model parameter.

Parameters:

sigma: float

The standard deviation of the distribution.

mu: float

The expectation of the distribution.

logprob(theta)

This method returns the log probability of the distribution for the given parameter theta.

Parameters:

theta: float

The parameter value for which the likelihood is to be computed.

Returns:

The value of the likelihood for the input parameter.

transform_prior(quantile)

This method returns the parameter value defined from the quantile of the cumulative distribution defined in the prior. It effectively converts an interval from 0-1 to one that samples the prior distribution and is used by nested sampling to define prior likelihoods.

Parameters:

quantile: float, 0-1

The quantile for which the prior is to be computed

Returns:

prior: float

The parameter value for the chosen parameter prior distribution.

class ndspec.SamplingUtils.priorLogNormal(mu, sigma)

This class is used to compute a lognormal prior distribution during Bayesian sampling, for a given model parameter.

Parameters:

sigma: float

The standard deviation of the distribution.

mu: float

The expectation of the distribution.

logprob(theta)

This method returns the log probability of the distribution for the given parameter theta.

Parameters:

theta: float

The parameter value for which the likelihood is to be computed.

Returns:

The value of the likelihood for the input parameter.

transform_prior(quantile)

This method returns the parameter value defined from the quantile of the cumulative distribution defined in the prior. It effectively converts an interval from 0-1 to one that samples the prior distribution and is used by nested sampling to define prior likelihoods.

Parameters:

quantile: float, 0-1

The quantile for which the prior is to be computed

Returns:

prior: float

The parameter value for the chosen parameter prior distribution.

class ndspec.SamplingUtils.priorTruncNormal(mu, sigma, min, max)

This class is used to compute a truncated, normal prior distribution during Bayesian sampling, for a given model parameter.

Parameters:

sigma: float

The standard deviation of the distribution.

mu: float

The expectation of the distribution.

min: float

The minimum value below which the distribution is truncated.

max: float

The maximum value above which the distribution is truncated.

logprob(theta)

This method returns the log probability of the distribution for the given parameter theta.

Parameters:

theta: float

The parameter value for which the likelihood is to be computed.

Returns:

The value of the likelihood for the input parameter.

transform_prior(quantile)

This method returns the parameter value defined from the quantile of the cumulative distribution defined in the prior. It effectively converts an interval from 0-1 to one that samples the prior distribution and is used by nested sampling to define prior likelihoods.

Parameters:

quantile: float, 0-1

The quantile for which the prior is to be computed

Returns:

prior: float

The parameter value for the chosen parameter prior distribution.

ndspec.SamplingUtils.nested_sampling_priors(quantile_cube)

This function samples the prior distribution for the parameters stored in the sampling_params global variable. This function is to be used to sample the prior distribution when using nested sampling algorithms. Given N parameters, the function maps an N-dimensional unitary cube to the corresponding quantile in the distribution for each prior. For example, given a parameter with a Gaussian prior centered at 2, and another parameter with a uniform prior from 3 to 5, the quantile cube (0.5,0.5) would return (2, 4) - the center of each distribution.

Parameters:

quantile_cube: np.array(float), 0-1

An array containing the quantile for each prior distribution from which to sample a new value

Returns:

params: np.array(float)

An array containing the parameter values sampled from the priors, given the input quantiles.

ndspec.SamplingUtils.log_priors(theta, prior_dict)

This function computes the total log-probability of a set of priors, given a st of input parameter values. This function is called automatically within the likelihood methods labelled “mcmc”.

Parameters:

theta: np.array(float)

An array of parameter values for which to compute the priors

prior_dict: dictionary

A dictionary of prior objects, each containing a method called .logprob which returns the log-probability given the input parameter value

Returns:

logprior: float

A float containing the log-probability of the set of parameters, given their priors.

ndspec.SamplingUtils.sampling_cash_likelihood(theta)

This function computes the log-likelihood of Poisson-distributed data excluding priors, for a given set of parameter values theta. This is the likelihood that should be passed to nested sampling algorithms, which evaluate the priors separately from the likelihood. It requires the global variables sampling_names, sampling_params, sampling_data, sampling_noise,

sampling_exp, and sampling_bins beforehand.

Parameters:

theta: np.array(float)

An array of parameter values for which to compute the log likelihood.

Returns:

likelihood: float

The value of the summed Cash log-likelihood for the given parameter values.

ndspec.SamplingUtils.mcmc_cash_likelihood(theta)

This function computes the log-likelihood of Poisson-distributed data, and including priors, for a given set of parameter values theta. This is the likelihood that should be passed to MCMC sampling algorithms, which evaluate the priors together with the likelihood. It requires the global variables sampling_priors, sampling_names, sampling_params, sampling_data, sampling_noise, sampling_exp, and sampling_bins beforehand.

Parameters:

theta: np.array(float)

An array of parameter values for which to compute the log likelihood.

Returns:

likelihood: float

The value of the summed Cash log-likelihood for the given parameter values.

ndspec.SamplingUtils.sampling_gaussian_likelihood(theta)

This function computes the log-likelihood, using a Gaussian distribution and including priors, for a given set of parameter values theta. This is the likelihood that should be passed to nested sampling algorithms, which evaluate the priors separately from the likelihood. It requires the user to have set the global variables sampling_names, sampling_params, sampling_data, sampling_data_err and sampling_model beforehand.

Parameters:

theta: np.array(float)

An array of parameter values for which to compute the log likelihood.

Returns:

likelihood: float

The value of the chi-square log-likelihood for the given parameter values.

ndspec.SamplingUtils.mcmc_gaussian_likelihood(theta)

This function computes the log-likelihood, using a Gaussian distribution and including priors, for a given set of parameter values theta. It requires the user to have set the global variables sampling_priors, sampling_names, sampling_params, sampling_data, sampling_data_err and sampling_model beforehand.

Parameters:

theta: np.array(float)

An array of parameter values for which to compute the log likelihood.

Returns:

likelihood: float

The value of the chi-square log-likelihood for the given parameter values.

ndspec.SamplingUtils.process_emcee(sampler, labels=None, discard=2000, thin=100, values=None, get_autocorr=True)

Given a sampler emcee EnsamleSampler object, this function calculates and prints the autocorrelation length, and plots the trace plots of the walkers, the acceptance fraction, and the corner plot for the posteriors.

This function is meant for a quick look at the output of a chain, rather than for publication quality plots. All the plots produced by this function have more customization options than the default ones used here.

Parameters:

sampler: emcee.EnsamleSampler

The sampler from which to plot the data

labels: list(str)

A list of strings to use for naming the parameters in both the trace and corner plots

discard: int, default 2000

The number of steps used to define the burn-in period

thin: int, default 100

Use one every “thin” steps in the chain. Used to make plots clearer.

values: np.array(float), default None

An array of parameter values used to show the best fit or “true” value of each parameter in the corner plot.

Xspec library

class ndspec.XspecInterface.ModelInterface(lib_path, pars_path)

This class allows users to load a library file containing Xspec-compatible models (including the entire library that comes with a typical HEASOFT installation), initialize them into the class objects as class methods, and evaluate them in their own Python code. ModelInterface serves as the parent for two more classes which can interface with models in Fortran (including the default Xspec library) or C respectively.

Attributes:

models_info: dict

A dictionary to store model information. The keywords store the initialized model names (models_info[‘nthcomp’]), the model type and parameters (e.g. models_info[[‘nthcomp’][‘type’]), and the names, minimum and maximum parameter values, and units of each (e.g. models_info[[‘nthcomp’][‘parameters’][‘kTe’][‘unit’] = “keV”).

lib: DLL ctype

The compiled .so (if using a Linux system) or .dylib (if using MacOS) library file loaded. By default, the code looks for the Xspec models at path_to_heasoft/Xspec/platform_name/lib/libXSFunctions.so (or dylib) file produced by an Xspec installation, and it contains (among other functions) the entire library of Xspec spectral models.

_all_info: dict

A dictionary containing the information of every model in the loaded library, regardless of whether the user has initialized it for use or not. The structure is identical to models_info.

check_param_values(model_name, params)

This method ensures that the input parameters for a given model evaluation do not exceed its allowed bounds, by comparing the input values with the minimum and maximum stored in the models_info dictionary

Parameters:

model_name: str

A string containing the name of the model being computed

params: np.array, dtype=float32

An array of parameter values stored as float32 being used in the model computation

Output:

test_pars: bool

A boolean which returns false if one of the number of parameters being passed is incorrect, or if one of the values is out of bounds. In the latter case, the model evaluation returns NaN. If none of the above happen, the method evalutes to True.

load_models(models)

This method allows users to initialized multiple models simultaneously by passing a dictionary with model names and calling functions.

Parameters:

models: dict

A dictionary whose keys are identical to the function names users want to intialize in the class.

parse_models(input_file)

This method parses the Xspec file with the model name, type, parameters values and units, etc, and stores the necessary information in the _all_info dictionary.

Parameters:

input_file: str

A path to the model file to be parsed. By default, the method looks for the Xspec model list in path_to_heasoft/Xspec/src/manager/.

Output:

models_info: dict

A dictionary containing model names and types, parameter values, minimum nad maximum bounds, and parameter units for all models in the library.

print_model_info()

This method prints to terminal a list of all the models that are currently initialized and ready for use, as well as their model type, parameter names, as well as default/min/max values and units.

class ndspec.XspecInterface.FortranInterface(lib_path=None, pars_path=None)

This class implements the methods inherited from ModelInterface to import and execute Xspec-compatible Fortran models. Users can use it either for their own code, or to load the entire Xspec library that is included in a typical HEASOFT installation.

add_model(func, symbol=None)

This method initializes a given model by adding it to the library object as one of its methods - for example:

def powerlaw(ear, params):

pass

lib.add_model(powerlaw) model = lib.powerlaw(arguments)

Parameters:

func: function

An empty function with the same name and input parameters as the model to be added to the library object.

symbol: str, optional

A string containing the name of the function in the library to be called. If not provided, it defaults to the function name. Generally, a user does not need to provide this argument, as the function name is usually sufficient to identify the model in the library. However, if a user wants to use a different name for the function they are calling via ndspec, they can provide the original function name here.

initialize_heasoft()

This method calls the fninit HEASOFT function, which initializes cross sections and abundances and is required to correctly evaluate Xspec models outside of the Xspec command interface.

class ndspec.XspecInterface.CInterface(lib_path, pars_path)

This class implements the methods inherited from ModelInterface to import and execute Xspec-compatible C models. It is meant to be used with custom models not immediately available with Xspec installations, such as Relxill and its various flavours.

add_model(func)

This method initializes a given model by adding it to the library object as one of its methods - for example:

def powerlaw(ear, params):

pass

lib.add_model(powerlaw) model = lib.powerlaw(arguments)

Parameters:

func: function

An empty function with the same name and input parameters as the model to be added to the library object.

Model library

ndspec.Models.lorentz(array, params, grid_edges=False)

This model is a Lorentzian function, defined identically to Uttley and Malzac 2023. The input parameters are:

array: the array over which the Lorentzian is to be computed

f_pk: the peak frequency of the Lorentzian

q: the q-factor of the Lorentzian

rms: the normalization of the Lorentzian

grid_edges: specifies whether the input array contains all the edges of a binned grid (identically to xspec), or the grid midpoints

ndspec.Models.cross_lorentz(array1, array2, params, grid_edges=False)

This model is a complex Lorentzian function, defined identically to Uttley and Malzac 2023, and shifted by a fixed phase, defined identically to Mendez et al. 2023. The input parameters are:

array: the array over which the Lorentzian is to be computed

f_pk: the peak frequency of the Lorentzian

q: the q-factor of the Lorentzian

rms: the normalization of the Lorentzian

phase: the phase lag associated with the Lorentzian

ndspec.Models.powerlaw(array, params, grid_edges=False)

This model is a standard power-law. The input parameters are:

array: the array grid over which to compute the power-law

norm: the normalization of the powerlaw

slope: the slope over the powerlaw. Unlike in Xspec, this parameter does not implicitely assume a minus sign; it must be specified by the user.

ndspec.Models.brokenpower(array, params, grid_edges=False)

This model is a smoothly broken powerlaw, defined identically to eq. 10 in Ghisellini and Tavecchio 2009. The input parameters are:

array: the array over which to compute the broken powerlaw

norm: the normalization of the broken powerlaw

slope1: the slope of the broken powerlaw before the break

slope2: the slope of the broken powerlaw after the break

brk: the location of the break in the powerlaw

ndspec.Models.gaussian(array, params, grid_edges=False)

This model is a Gaussian function. The input parameters are:

array: the array over which the Gaussian is defined

center: the centroid of the Gaussian

width: the width of the Gaussian

gauss_norm: the normalization of the Gaussian

ndspec.Models.bbody(array, params, grid_edges=False)

This model is a constant black body. The input parameters are:

array: the array over which the spectrum is defined

norm: the normalization of the black body, defined identically to that of the Xspec model

temp: the temperature in keV

ndspec.Models.varbbody(array, params, grid_edges=False)

This model is a variable black body, defined identically to Uttley and Malzac 2023. The input parameters are:

array: the array over which the spectrum is defined

norm: the normalization of the black body, defined identically to that of the Xspec model

temp: the temperature in keV

ndspec.Models.gauss_fred(array1, array2, params, return_full=False, grid_edges=False)

This is a two-dimensional model for an impulse response function. The time dependence is a fast rise, exponential decay pulse. The dependence over the second axis (typically energy) is a Gaussian line narrowing over time following a powerlaw. The total model is the product of the two dependences. The input paramters are:

array1: the time over which the pulse is defined

array2: the second direction over which the model is defined

norm: the total model normalization

width: the initial width of the Gaussian

center: the centroid of the Gaussian

rise_t: the rise pulse timescale

decay_t: the decay pulse timescale

decay_w: the slope of the energy width powerlaw decay

return_full: a boolean to choose whether to return just the 2d model (done by default), or the additional projections over the two model axis

ndspec.Models.gauss_bkn(array1, array2, params, return_full=False)

This is a two-dimensional model for an impulse response function. The time dependence is a smoothly broken powerlaw pulse. The dependence over the second axis (typically energy) is a Gaussian line narrowing over time following a powerlaw. The total model is the product of the two dependences. The input paramters are:

array1: the time over which the pulse is defined

array2: the second direction over which the model is defined

norm: the total model normalization

center: the centroid of the Gaussian

width: the initial width of the Gaussian

rise_slope: the rise pulse slope

decay_slope: the decay pulse slope

break_time: the time at which the broken powerlaw changes from rise to decay slope

decay_w: the slope of the energy width powerlaw decay

return_full: a boolean to choose whether to return just the 2d model (done by default), or the additional projections over the two model axis

ndspec.Models.bbody_fred(array1, array2, params, return_full=False)

This is a two-dimensional model for an impulse response function. The time dependence is a fast rise, exponential decay pulse. The dependence over the second energy is a variable black body, cooling over time following a powerlaw. The total model is the product of the two dependences. The input paramters are:

array1: the time over which the pulse is defined

array2: the second direction over which the model is defined

norm: the total model normalization

temp: the initial temperature

rise_slope: the rise pulse slope

decay_slope: the decay pulse slope

break_time: the time at which the broken powerlaw changes from rise to decay slope

decay_temp: the slope of the temperature powerlaw decay

return_full: a boolean to choose whether to return just the 2d model (done by default), or the additional projections over the two model axis

ndspec.Models.bbody_bkn(array1, array2, params, return_full=False)

This is a two-dimensional model for an impulse response function. The time dependence is a smoothly broken powerlaw pulse. The dependence over the second energy is a variable black body, cooling over time following a powerlaw. The total model is the product of the two dependences. The input paramters are:

array1: the time over which the pulse is defined

array2: the second grid over which the model is defined

norm: the total model normalization

temp: the initial temperature

rise_slope: the rise pulse slope

decay_slope: the decay pulse slope

break_time: the time at which the broken powerlaw changes from rise to decay slope

decay_temp: the slope of the temperature powerlaw decay

return_full: a boolean to choose whether to return just the 2d model (done by default), or the additional projections over the two model axis

ndspec.Models.pivoting_pl(array1, array2, params)

This is a pivoting power-law model similar to that implemented in reltrans (Mastroserio et al. 2021). The main difference is that this implementation expresses the dependence of the paramters gamma and phi_ab (in the paper above) explicitely. The input paramters are:

array1: the Fourier frequencies over which to compute the model

array2: the second direction (typically energy) over which to compute the model

norm: the model normalzation

pl_index: the slope over the powerlaw

gamma_0: the gamma parameter in Mastroserio et al. 2021, defined at a frequency nu_0

gamma_slope: the dependence of the gamma parameter with Fourier frequency, which is assumed to be log-linear

phi_0: the phi_AB parameter in Mastroserio et al. 2021, defined at a frequency nu_0

nu_0: the initial frequency from which the pivoting parameters are defined

ndspec.Models.plot_2d(xaxis, yaxis, impulse_2d, impulse_x, impulse_y, xlim=[0.0, 400.0], ylim=[0.1, 10.5], xlog=False, ylog=False, return_plot=False, normalize_en=True)

A simple automated plotter for the impulse response function models above. The input parameters are:

xaxis, yaxis: the two grids over which the model is defined

impulsed_2d: the two-dimensional model to plot

impulse_x,impulse_y: the projections of the model over the x/y axis

xlim,ylim: the limits of the x/y axis to show in the plot

xlog,ylog: booleans to switch between linera and log scales in each axis

return_plot: boolean to return the figure object for storage/saving

normalize_en: boolean to multiply the energy dependence (on the y axis) by the y axis values squared. Useful to highlight the model energy dependence.

Simulator utilities

ndspec.Simulator.simulate_lightcurve(psd_obj, obs_time, dt, countrate, rms=None, params=None)

This function is used to simulate a lightcurve of the set model for a given set of parameters for a given timespan at a given time resolution. By default, this will simulate using the specified model and model parameters and will self-consistently calculate the root mean square of the lightcurve from the power spectrum model. The count rate must always be provided, as the mean flux of the lightcurve cannot be determined from the model alone. This method does not take into account the instrument response, and therefore all lightcurves simulated with this method will be mono-energetic and will not include any energy-dependent effects.

Practically, this method is a wrapper for the stingray simulator module (https://docs.stingray.science/en/stable/api.html#simulator).

Parameters:

psd_obj: ndspec.PowerSpectrum

An instance of the PowerSpectrum class, which contains the model to use for simulating the lightcurve. This object is used to evaluate the power spectrum at the frequencies measurable by the observation time and time resolution.

obs_time: float

The total observation time, in seconds, over which the lightcurve is simulated. This is used to determine the number of bins in the simulated lightcurve.

dt: float

The time resolution of the simulated lightcurve, in seconds. This determines the time binning of the simulated lightcurve.

countrate: float

The mean count rate of the simulated lightcurve, in counts per second. This is used to set the mean flux of the simulated lightcurve.

rms: float, default None

The root mean square of the simulated lightcurve, which is used to set the variability of the simulated lightcurve. By default, the rms is calculated from the power spectrum model. If a specific rms value is provided, it will be used instead.

params: lmfit.Parameters, default None

The parameter values to use in evaluating the model. If none are provided, the model_params attribute of the FitTimeAvgSpectrum is used.

Returns:

lightcurve: stingray.lightcurve.Lightcurve

The resulting lightcurve object, containing the simulated lightcurve of the model evaluated over the given Fourier frequency array, for the given input parameters.

ndspec.Simulator.simulate_time_lags(fitobj, ref_Elo, ref_Ehi, sub_Elo, sub_Ehi, Texp, coh2, pow, time_avg_model, bkg_file_path=None, params=None)

This method will simulate a lag spectrum based on the model. The model must be able to evaluate the cross spectrum and a background file must be provided. You must also provide the energy bounds of the reference and subject bands, the exposure time, the coherence squared value, and the power in rms/mean$^2$/Hz units ($lpha_{

u}$), as well as a time-averaged model that shares the same physical

assumptions (and thus share relevant parameters) as the set cross-spectrum model. The method will return a one-dimensional array containing the simulated lags for each energy channel, based on the model and the specified parameters.

For further explanation for how the lag spectrum is simulated, see section 3 of Ingram et al. 2022, https://ui.adsabs.harvard.edu/abs/2022MNRAS.509..619I/abstract.

fitobj: ndspec.FitCrossSpectrum

An instance of the FitCrossSpectrum class, which contains the model to use for simulating the lag spectrum. The FitCrossSpectrum object must have a model defined, a frequency grid set, and an instrument response set before calling this function.

ref_Elo: float

The lower energy bound of the reference band in keV.

ref_Ehi: float

The upper energy bound of the reference band in keV.

sub_Elo: float

The lower energy bound of the subject band in keV.

sub_Ehi: float

The upper energy bound of the subject band in keV.

Texp: float

The exposure time in seconds for which the lag spectrum is simulated.

coh2: float

The coherence squared value, which is a measure of the correlation between the reference and subject bands.

pow: float

The power in rms/mean$^2$/Hz units ($lpha_{

u}$), which is used to scale the

cross spectrum.

time_avg_model: lmfit.Model or lmfit.CompositeModel

A model that evaluates the time-averaged power spectrum, which is used to calculate the background noise in the reference and subject bands. This model should share the same physical assumptions as the model used to evaluate the cross spectrum and share all relevant physical parameters

bkg_file_path: str, optional

The path to the background file containing the background counts for each energy channel. If not provided, the method will default to the background file already set in the class. A background file must be provided to simulate the lag spectrum.

params: lmfit.Parameters, optional

The parameters to use for evaluating the model. If not provided, the default parameters stored in the model_params attribute will be used.

lagsim: np.array(float)

A one-dimensional array containing the simulated lags for each energy channel, based on the model and the specified parameters. The size of this array is equal to the number of energy channels defined in the instrument response matrix.

ndspec.Simulator.simulate_time_averaged(res_obj, model, params, exposure_time=None)

This method simulates a time-averaged spectrum given a set of parameters, by evaluating the model and folding it through the response. It is used to generate synthetic spectra for testing purposes.

Parameters:

res_obj: ndspec.ResponseMatrix

An instance of the ResponseMatrix class, which contains the response matrix of the instrument for which to simulate the data

model: lmfit.model or lmfit.compositemodel

An instance of an LMFit model object which contains the model to use to simulate the spectrum

params: lmfit.Parameters, default None

The parameter values to use in evaluating the model. If none are provided, the model_params attribute is used.

exposure_time: float, default None

The exposure time to use for the simulation. If None, the exposure time stored in the response matrix is used. This is used to convert the model count rate to expected counts in each channel.

ear: np.array(float), default None

The array of photon energy channel edges over which to evaluate the model. If none are provided, the same grid contained in the instrument response is used.

Returns:

simulated_spectrum: np.array(float)

The simulated spectrum evaluated over the noticed energy channels and Poisson sampled. The spectrum is in units of counts/channel.