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)

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)

The quantity from which to calculate the power spectrum.

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

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

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: None

Work in progress; currently the software defaults to chi squared 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.

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 “delchi”, it returns the contribution of each energy channel to the total chi squared.

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_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.

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 Aadjusts 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.

FitPowerSpectrum Class

class ndspec.FitPowerSpectrum.FitPowerSpectrum

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: None

Work in progress; currently the software defaults to chi squared 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='delchi', 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=”delchi”

The units to use for the residuals. If residuals=”delchi”, the plot shows the residuals in units of data-model/error; if residuals=”ratio”, the plot instead uses units of data/model.

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

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: None

Work in progress; currently the software defaults to chi squared 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.

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.

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.

eval_model(params=None, energ=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.

energ: np.array(float), default None

The the photon energies over which to evaluted 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', 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.

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='data', residuals='delchi', 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.

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=”delchi”

The units to use for the residuals. If residuals=”delchi”, the plot shows the residuals in units of data-model/error; if residuals=”ratio”, the plot instead uses units of data/model.

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)

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).

FitCrossSpectrum Class

class ndspec.FitCrossSpectrum.FitCrossSpectrum

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: None

Work in progress; currently the software defaults to chi squared 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.

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).

_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.

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, use_phase=False, residuals='delchi', 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=”delchi”

The units to use for the residuals. If residuals=”delchi”, the plot shows the residuals in units of data-model/error; if residuals=”ratio”, the plot instead uses units of data/model.

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='delchi', 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=”delchi”

The units to use for the residuals. If residuals=”delchi”, the plot shows the residuals in units of data-model/error; if residuals=”ratio”, the plot instead uses units of data/model.

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(value)

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:

value: 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.

renorm_phases(value)

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:

value: 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_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.

Emcee sampling functions

ndspec.EmceeUtils.set_emcee_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 emcee_priors; therefore, users should never re-use the variable name emcee_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.EmceeUtils.set_emcee_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 emcee_model; therefore, users should never re-use the variable name emcee_model in their code.

Parameters:

fitobj: ndspec.Fit…Object or ndspec.JointFit

Object containing the data, specified model, and parameters

ndspec.EmceeUtils.set_emcee_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 emcee_data and emcee_data_err; therefore, users should never re-use the variable names emcee_data and emcee_data_err in their code.

Parameters:

fitobj: ndspec.Fit…Object or ndspec.JointFit

Object containing the data, specified model, and parameters

ndspec.EmceeUtils.set_emcee_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 emcee_names, emcee_values and emcee_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.

class ndspec.EmceeUtils.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.

class ndspec.EmceeUtils.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 log(x) is uniform, then p(log(x)) = const, p(x) = p(log(x))*dlog(x)/dx = const/x. Therefore, the log-probability is (minus a constant) log(p(x)) = log(1/x) = -log(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.

class ndspec.EmceeUtils.priorNormal(sigma, mu)

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.

class ndspec.EmceeUtils.priorLogNormal(sigma, mu)

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.

ndspec.EmceeUtils.log_priors(theta, prior_dict)

This function computes the total log-probability of a set of priors, given a st of input parameter values.

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.EmceeUtils.chi_square_likelihood(theta)

This function computes the log-likelihood, using the chi-square statistic and including priors, for a given set of parameter values theta. It requires the user to have set the global variables emcee_priors, emcee_names, emcee_params, emcee_data, emcee_data_err and emcee_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.EmceeUtils.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.

Model library

ndspec.models.lorentz(array, params)

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

ndspec.models.cross_lorentz(array1, array2, params)

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)

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)

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)

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)

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)

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)

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.