py4py.reverb

Contains the class used to create and manipulate reverberation maps from Python output files.

Example

For an existing delay output file called ‘qso.delay_dump’, to generate a TF plot for the C4 line for a specific spectrum, with axis of velocity offset vs days, you would do:

qso_conn = open_database('qso')
tf_c4_1 = TransferFunction(
    qso_conn, continuum=1e43, wave_bins=100, delay_bins=100, filename='qso_c4_spectrum_1'
)
tf_c4_1.spectrum(1).line(443).run()
tf_c4_1.plot(velocity=True, days=True)

Given database queries can take a long time, it is advisable to pickle a TF that has been run so you can access it later on. Note, however: Once a TF has been restored from a pickle, you can no longer change the filters and re-run:

with open('qso_c4_spectrum_1', 'wb') as file:
    pickle.dump(tf_c4_1, file)

py4py.reverb.timeseries

py4py.reverb.output

`SECONDS_PER_DAY`	Constant used for rescaling data, that is probably superfluous and already present in Astropy
`Base`	Base class declared dynamically to bind to SQLalchemy
`kep_sey`	The default Keplerian outline settings for the Seyfert model
`kep_qso`	The default Keplerian outline settings for the QSO model

`TransferFunction`	Used to create, store and query emissivity and response functions
`Spectrum`	The SQLalchemy table for the spectra. Unused.
`Origin`	The SQLalchemy table for the photon origins. Unused.
`Photon`	SQLalchemy class for a photon. Why are all the properties capitalised?

`calculate_fwhm`(→ float)	Calculate FWHM from arrays
`calculate_midpoints`(→ numpy.typing.NDArray[numpy.floating])	Converts bin boundaries into midpoints
`keplerian_velocity`(→ float)	Calculates Keplerian velocity at given radius
`doppler_shift_wave`(→ float)	Converts passed line and velocity into red/blue-shifted wavelength
`doppler_shift_vel`(→ float)	Converts passed red/blue-shifted wave into velocity
`calculate_delay`(→ float)	Delay relative to continuum for emission from a point on the disk.
`open_database`(→ sqlalchemy.engine.Engine)	Open or create a SQL database

Package Contents

py4py.reverb.calculate_fwhm(midpoints: numpy.typing.NDArray[numpy.floating], vals: numpy.typing.NDArray[numpy.floating]) → float

Calculate FWHM from arrays

Taken from http://stackoverflow.com/questions/10582795/finding-the-full-width-half-maximum-of-a-peak I don’t think this can cope with being passed a doublet or an array with no peak within it. Doublets will calculate FWHM from the HM of both!

Parameters:

midpoints (numpy.typing.NDArray[numpy.floating]) – Array of bin midpoints
vals (numpy.typing.NDArray[numpy.floating]) – Array of bin values

Returns:

FWHM of the peak (should it exist!)

Return type:

float

py4py.reverb.calculate_midpoints(bins: numpy.typing.NDArray[numpy.floating]) → numpy.typing.NDArray[numpy.floating]

Converts bin boundaries into midpoints

Parameters:: bins (numpy.typing.NDArray[numpy.floating]) – Array of bin boundaries
Returns:: Array of bin midpoints (1 shorter!)
Return type:: numpy.typing.NDArray[numpy.floating]

py4py.reverb.keplerian_velocity(mass: float, radius: float) → float

Calculates Keplerian velocity at given radius

Parameters:

mass (float) – Object mass in kg
radius (float) – Orbital radius in m

Returns:

Orbital velocity in m/s

Return type:

float

py4py.reverb.doppler_shift_wave(line: float, vel: float) → float

Converts passed line and velocity into red/blue-shifted wavelength

Parameters:

line (float) – Line wavelength (any length unit)
vel (float) – Doppler shift velocity (m/s)

Returns:

Doppler shifted line wavelength (as above)

Return type:

float

py4py.reverb.doppler_shift_vel(line: float, wave: float) → float

Converts passed red/blue-shifted wave into velocity

Parameters:

line (float) – Base line wavelength (any length unit)
wave (float) – Doppler shifted line wavelength (as above)

Returns:

Speed of Doppler shift

Return type:

float

py4py.reverb.SECONDS_PER_DAY = 86400: Constant used for rescaling data, that is probably superfluous and already present in Astropy

py4py.reverb.calculate_delay(angle: float, phase: float, radius: float, days: bool = True) → float

Delay relative to continuum for emission from a point on the disk.

Calculate delay for emission from a point on a keplerian disk, defined by its radius and disk angle, to an observer at a specified angle.

Draw plane at r_rad_min out. Find x projection of disk position. Calculate distance travelled to that plane from the current disk position Delay relative to continuum is thus (distance from centre to plane) + distance from centre to point

Parameters:

angle (float) – Observer angle to disk normal, in radians
phase (float) – Rotational angle of point on disk, in radians. 0 = in line to observer
radius (float) – Radius of the point on the disk, in m
days (bool) – Whether the timescale should be seconds or days

Returns:

Delay relative to continuum

Return type:

float

class py4py.reverb.TransferFunction(database: sqlalchemy.engine.Connection, filename: str, continuum: float, wave_bins: int = None, delay_bins: int = None, template: TransferFunction = None, template_different_line: bool = False, template_different_spectrum: bool = False)

Used to create, store and query emissivity and response functions

Initialises the TF, optionally by templating off another TF.

Sets up all the basic properties of the TF that are required to create it. It must be .run() to query the DB before it can itself be queried. If templating, it applies all the same filters that were applied to the template TF, unless explicitly told not to. Filters don’t overwrite! They stack. So you can’t simply call .line() to change the line the TF corresponds to if its template was a different line, unless you specify that the template was of a different line.

Parameters:

database (sqlalchemy.engine.Connection) – The database to be queried for this TF.
filename (string) – The root filename for plots created for this TF.
continuum (float) – The continuum value associated with this TF. Central source + disk luminosity.
wave_bins (int) – Number of wavelength/velocity bins.
delay_bins (int) – Number of delay time bins.
template (TransferFunction) – Other TF to copy all filter settings from. Will match delay, wave and velocity bins exactly.
template_different_line (bool) – Is this TF going to share delay & velocity bins but have different wavelength bins?
template_different_spectrum (bool) – Is this TF going to share all specified bins but be taken on photons from a different observer.

Todo

Consider making it impossible to apply filters after calling run().

__getstate__() → dict

Removes invalid data before saving to disk.

Returns:

Updated internal dict, with references to external,: session-specific database things, removed.

Return type:

dict

__setstate__(state: dict)

Restores the data from disk, and sets a flag to show this is a frozen TF.

Parameters:: state (dict) – The unpickled object dict..

_database

_session

_query

_delay_dynamic_range = None

_velocity = None

_line_list = None

_line_wave = None

_line_num = None

_delay_range = None

_continuum

_filename

_bins_wave_count = None

_bins_delay_count = None

_bins_vel = None

_bins_wave = None

_bins_delay = None

_emissivity = None

_response = None

_count = None

_wave_range = None

_spectrum = None

_unpickled = False

spectrum(number: int) → TransferFunction

Constrain the TF to photons from a specific observer

Parameters:: number (int) – Observer number from Python run
Returns:: Self, so filters can be stacked
Return type:: TransferFunction

line(number: int, wavelength: float) → TransferFunction

Constrain the TF to only photons last interacting with a given line

This includes being emitted in the specified line, or scattered off it

Parameters:

number (int) – Python line number. Will vary based on data file!
wavelength (float) – Wavelength of the line in angstroms

Returns:

Self, so filters can be stacked

Return type:

TransferFunction

velocities(velocity: float) → TransferFunction

Constrain the TF to only photons with a range of Doppler shifts

Parameters:: velocity (float) – Maximum doppler shift velocity in m/s. Applies to both positive and negative Doppler shift
Returns:: Self, so filters can be stacked
Return type:: TransferFunction

wavelengths(wave_min: float, wave_max: float) → TransferFunction

Constrain the TF to only photons with a range of wavelengths

Parameters:

wave_min (float) – Minimum wavelength in angstroms
wave_max (float) – Maximum wavelength in angstroms

Returns:

Self, so filters can be stacked

Return type:

TransferFunction

wavelength_bins(wave_range: numpy.ndarray) → TransferFunction

Constrain the TF to only photons with a range of wavelengths, and to a specific set of bins

Parameters:: wave_range (np.ndarray) – Array of bins to use
Returns:: Self, so filters can be stacked
Return type:: TransferFunction

lines(line_list: List[int]) → TransferFunction

Constrain the TF to only photons with a specific internal line number. This list number will be specific to the python atomic data file!

Parameters:: line_list (List[int]) – List of lines
Returns:: Self, so filters can be stacked
Return type:: TransferFunction

delays(delay_min: float, delay_max: float, days: bool = True) → TransferFunction

The delay range that should be considered when producing the TF.

Parameters:

delay_min (float) – Minimum delay time (in seconds or days)
delay_max (float) – Maximum delay time (in seconds or days)
days (bool) – Whether or not the delay range has been provided in days

Returns:

Self, so filters can be stacked

Return type:

TransferFunction

delay_dynamic_range(delay_dynamic_range: float) → TransferFunction

If set, the TF will generate delay bins to cover this dynamic range of responses, i.e. (1 - 10^-ddr) of the delays. So a ddr of 1 will generate photons with delays up to 1 - (1/10) = the 90th percentile of delays. ddr=2 will give up to the 99th percentile, 3=99.9th percentile, etc.

Arguably this is a bit of an ambiguous name

Parameters:: delay_dynamic_range (float) – The dynamic range to be used when
Returns:: Self, so filters can be stacked
Return type:: TransferFunction

cont_scatters(scat_min: int, scat_max: int | None = None) → TransferFunction

Constrain the TF to only photons that have scattered min-max times via a continuum scattering process (e.g. electron scattering).

Parameters:

scat_min (int) – Minimum number of continuum scatters
scat_max (Optional[int]) – Maximum number of continuum scatters, if desired

Returns:

Self, so filters can be stacked

Return type:

TransferFunction

res_scatters(scat_min: int, scat_max: int | None = None) → TransferFunction

Constrain the TF to only photons that have scattered min-max times via a resonant scattering process (e.g. line scattering).

Parameters:

scat_min (int) – Minimum number of resonant scatters
scat_max (Optional[int]) – Maximum number of resonant scatters, if desired

Returns:

Self, so filters can be stacked

Return type:

TransferFunction

filter(*args) → TransferFunction

Apply a SQLalchemy filter directly to the content.

Parameters:: args – The list of filter arguments
Returns:: Self, so filters can be stacked
Return type:: TransferFunction

response_map_by_tf(low_state: TransferFunction, high_state: TransferFunction, cf_low: float = 1.0, cf_high: float = 1.0) → TransferFunction

Creates a response function for this transfer function by subtracting two transfer functions bracketing it. Requires two other completed transfer functions, bracketing this one in luminosity, all with matching wavelength/velocity and delay bins.

Correction factors are there to account for things like runs that have been terminated early, e.g. if you request 100 spectrum cycles and stop (or Python dies) after 80, the total photon luminosity will only be 80/100. A correction factor allows you to bump this up. Arguably correction factors should be applied during the ‘run()’ method.

Parameters:

low_state (TransferFunction) – A full, processed transfer function for a lower-luminosity system.
high_state (TransferFunction) – A full, processed transfer function for a higher-luminosity system.
cf_low (float) – Correction factor for low state. Multiplier to the whole transfer function.
cf_high (float) – Correction factor for high state. Multiplier to the whole transfer function.

Returns:

Self, so plotting can be chained on.

Return type:

TransferFunction

fwhm(response: bool = False, velocity: bool = True)

Calculates the full width half maximum of the delay-summed transfer function, roughly analogous to the line profile. Possibly meaningless for the response function?

Parameters:

response (bool) – Whether to calculate the FWHM of the transfer or response function
velocity (bool) – Whether to return the FWHM in wavelength or velocity-space

Returns:

Full width at half maximum for the function.: If the function is a doublet, this will not work properly.

Return type:

float

Todo

Catch doublets.

delay(response: bool = False, threshold: float = 0, bounds: float = None, days: bool = False) → float | Tuple[float, float, float]

Calculates the centroid delay for the current data

Parameters:

response (bool) – Whether or not to calculate the delay from the response
threshold (float) – Exclude all bins with value < the threshold fraction of the peak value. Standard value used in the reverb papers was 0.8.
bounds (float) – Return the fractional bounds (i.e. bounds=0.25, the function will return [0.5, 0.25, 0.75]). Not implemented.
days (bool) – Whether to return the delay in days or seconds

Returns:

Centroid delay, and lower and upper fractional bounds if bounds keyword provided

Return type:

Union[float, Tuple[float, float, float]]

Todo

Implement fractional bounds. Should just be able to call the centroid_delay function!

delay_peak(response: bool = False, days: bool = False) → float

Calculates the peak delay for the transfer or response function, i.e. the delay at which the response is strongest.

Parameters:

response (bool) – Whether or not to calculate the peak transfer or response function.
days (bool) – Whether to return the value in seconds or days.

Returns:

The peak delay.

Return type:

float

run(scaling_factor: float = 1.0, limit: int = None, verbose: bool = False) → TransferFunction

Performs a query on the photon DB and bins it.

A TF must be run after all filters are applied and before any attempts to retrieve or process data from it. This can be a time-consuming call, on the order of 1 minute per GB of input file.

Parameters:

scaling_factor (float) – 1/Number of cycles in the spectra file
limit (int) – Number of photons to limit the TF to, for testing. Recommend testing filters on a small number of photons to begin with.
verbose (bool) – Whether to output exactly what the query is.

Returns:

Self, for chaining commands

Return type:

TransferFunction

Internal function used by response(), emissivity() and count()

Parameters:

array (np.ndarray) – Array to return value from
delay (Optional[float]) – Delay to return value for. Must provide this or delay_index.
delay_index (Optional[int]) – Delay index to return value for. Must provide this or delay.
wave (Optional[float]) – Wavelength to return value for

Returns:

Either a subset of the array if only delay is provided,: or the value of a single array element if delay and wavelength provided.

Return type:

Union[np.ndarray, float]

Todo

Allow for only wavelength to be provided?

response_total() → float

Returns the total response.

Returns:: Total response.
Return type:: float

delay_bins() → numpy.ndarray

Returns the range of delays covered by this TF.

Returns:: Array of the bin boundaries.
Return type:: np.ndarray

response(delay: float | None = None, wave: float | None = None, delay_index: int | None = None) → float | numpy.ndarray

Returns the responsivity in either one specific wavelength/delay bin, or all wavelength bins for a given delay.

Parameters:

delay (Optional[float]) – Delay to return value for. Must provide this or delay_index.
delay_index (Optional[int]) – Delay index to return value for. Must provide this or delay.
wave (Optional[float]) – Wavelength to return value for.

Returns:

Either the responsivity in one specific bin, or if wave is not specified: the counts in each wavelength bin at this delay

Return type:

Union[int, np.ndarray]

Todo

Allow for only wavelength to be provided?

emissivity(delay: float | None = None, wave: float | None = None, delay_index: int | None = None) → float | numpy.ndarray

Returns the emissivity in either one specific wavelength/delay bin, or all wavelength bins for a given delay.

Parameters:

delay (Optional[float]) – Delay to return value for. Must provide this or delay_index.
delay_index (Optional[int]) – Delay index to return value for. Must provide this or delay.
wave (Optional[float]) – Wavelength to return value for.

Returns:

Either the emissivity in one specific bin, or if wave is not specified: the counts in each wavelengthin bin at this delay

Return type:

Union[int, np.ndarray]

Todo

Allow for only wavelength to be provided?

count(delay: float | None = None, wave: float | None = None, delay_index: int | None = None) → int | numpy.ndarray

Returns the photon count in either one specific wavelength/delay bin, or all wavelength bins for a given delay.

Parameters:

delay (Optional[float]) – Delay to return value for. Must provide this or delay_index.
delay_index (Optional[int]) – Delay index to return value for. Must provide this or delay.
wave (Optional[float]) – Wavelength to return value for

Returns:

Either the count in one specific bin, or if wave is not specified: the counts in each wavelength bin at this delay

Return type:

Union[int, np.ndarray]

Todo

Allow for only wavelength to be provided?

transfer_function_1d(response: bool = False, days: bool = True) → numpy.ndarray

Collapses the 2-d transfer/response function into a 1-d response function, and returns the bin midpoints and values in each bin for plotting.

Parameters:

response (bool) – Whether or not to return the response function data
days (bool) – Whether the bin midpoints should be in seconds or days

Returns:

A [bins, 2]-d array containing the midpoints of the delay bins,: and the value of the 1-d transfer or response function in each bin.

Return type:

np.ndarray

plot(log: bool = False, normalised: bool = False, rescaled: bool = False, velocity: bool = False, name: str = None, days: bool = True, response_map=False, keplerian: dict = None, dynamic_range: int = None, rms: bool = False, show: bool = False, max_delay: float | None = None, format: str = '.eps', return_figure: bool = False) → TransferFunction | matplotlib.figure.Figure

Takes the data gathered by calling ‘run’ and outputs a plot

Parameters:

log (bool) – Whether the plot should be linear or logarithmic.
normalised (bool) – Whether or not to rescale the plot such that the total emissivity = 1.
rescaled (bool) – Whether or not to rescale the plot such that the maximum emissivity = 1.
velocity (bool) – Whether the plot X-axis should be velocity (true) or wavelength (false).
name (Optional[str]) – The file will be output to ‘tf_filename.eps’. May add the ‘name’ component to modify it to ‘tf_filename_name.eps’. Useful for adding e.g. ‘c4’ or ‘log’.
days (bool) – Whether the plot Y-axis should be in days (true) or seconds (false).
response_map (bool) – Whether to plot the transfer function map or the response function.
keplerian (Optional[dict]) – A dictionary describing the profile of a keplerian disk, the bounds of which will be overlaid on the plot. Arguments include angle (float) - Angle of disk to the observer, mass (float) - Mass of the central object in M_sol, radius (Tuple(float, float)) - Inner and outer disk radii, in $r_{ISCO}$. include_minimum_velocity - Whether or not to include the outer disk velocity profile (default no).
dynamic_range (Optional[int]) – If the plot is logarithmic, the dynamic range the colour bar should show. If not provided, will attempt to use the base dynamic range property, otherwise will default to showing 99.9% of all emissivity.
max_delay (Optional[float]) – The optional maximum delay to plot out to.
rms (bool) – Whether or not the line profile panel should show the root mean squared line profile.
show (bool) – Whether or not to display the plot to screen.
format (str) – The output file format. .eps by default.
return_figure (bool) – If true, return the figure instead of platting it.

Returns:

Self, for chaining outputs

Return type:

TransferFunction

py4py.reverb.open_database(file_root: str, user: str = None, password: str = None, batch_size: int = 25000) → sqlalchemy.engine.Engine

Open or create a SQL database

Will open a SQLite DB if one already exists, otherwise will create one from file. Note, though, that if the process is interrupted the code cannot intelligently resume - you must delete the half-written DB!

Parameters:

file_root (string) – Root of the filename (no ‘.db’ or ‘.delay_dump’)
user (string) – Username. Here in case I change to PostgreSQL
password (string) – Password. Here in case I change to PostgreSQL
batch_size (int) – Number of photons to stage before committing. If too low, file creation is slow. If too high, get out-of-memory errors.

Returns:

Connection to the database opened

Return type:

sqlalchemy.engine.Engine

py4py.reverb.Base: Base class declared dynamically to bind to SQLalchemy

class py4py.reverb.Spectrum

Bases: Base

The SQLalchemy table for the spectra. Unused. Could be removed but will break backward compatibility. Information required for this is not stored in the output files.

# Todo: Implement or remove this table.

__tablename__ = 'Spectra'

id

angle

class py4py.reverb.Origin

Bases: Base

The SQLalchemy table for the photon origins. Unused. Could be removed but will break backward compatibility. Information required for this is not stored in the output files.

# Todo: Implement or remove this table