py4py.reverb ============ .. py:module:: py4py.reverb .. autoapi-nested-parse:: Contains the class used to create and manipulate reverberation maps from Python output files. .. rubric:: 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) .. toctree:: :maxdepth: 1 /py_progs/py4py_auto/py4py/reverb/timeseries/index .. toctree:: :maxdepth: 1 /py_progs/py4py_auto/py4py/reverb/output/index .. autoapisummary:: py4py.reverb.SECONDS_PER_DAY py4py.reverb.Base py4py.reverb.kep_sey py4py.reverb.kep_qso .. autoapisummary:: py4py.reverb.TransferFunction py4py.reverb.Spectrum py4py.reverb.Origin py4py.reverb.Photon .. autoapisummary:: py4py.reverb.calculate_fwhm py4py.reverb.calculate_midpoints py4py.reverb.keplerian_velocity py4py.reverb.doppler_shift_wave py4py.reverb.doppler_shift_vel py4py.reverb.calculate_delay py4py.reverb.open_database Package Contents ---------------- .. py:function:: 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! :param midpoints: Array of bin midpoints :param vals: Array of bin values :returns: FWHM of the peak (should it exist!) .. py:function:: calculate_midpoints(bins: numpy.typing.NDArray[numpy.floating]) -> numpy.typing.NDArray[numpy.floating] Converts bin boundaries into midpoints :param bins: Array of bin boundaries :returns: Array of bin midpoints (1 shorter!) .. py:function:: keplerian_velocity(mass: float, radius: float) -> float Calculates Keplerian velocity at given radius :param mass: Object mass in kg :type mass: float :param radius: Orbital radius in m :type radius: float :returns: Orbital velocity in m/s :rtype: float .. py:function:: doppler_shift_wave(line: float, vel: float) -> float Converts passed line and velocity into red/blue-shifted wavelength :param line: Line wavelength (any length unit) :type line: float :param vel: Doppler shift velocity (m/s) :type vel: float :returns: Doppler shifted line wavelength (as above) :rtype: float .. py:function:: doppler_shift_vel(line: float, wave: float) -> float Converts passed red/blue-shifted wave into velocity :param line: Base line wavelength (any length unit) :type line: float :param wave: Doppler shifted line wavelength (as above) :type wave: float :returns: Speed of Doppler shift :rtype: float .. py:data:: SECONDS_PER_DAY :value: 86400 Constant used for rescaling data, that is probably superfluous and already present in Astropy .. py:function:: 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 :param angle: Observer angle to disk normal, in radians :type angle: float :param phase: Rotational angle of point on disk, in radians. 0 = in line to observer :type phase: float :param radius: Radius of the point on the disk, in m :type radius: float :param days: Whether the timescale should be seconds or days :type days: bool :returns: Delay relative to continuum :rtype: float .. py:class:: 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. :param database: The database to be queried for this TF. :type database: sqlalchemy.engine.Connection :param filename: The root filename for plots created for this TF. :type filename: string :param continuum: The continuum value associated with this TF. Central source + disk luminosity. :type continuum: float :param wave_bins: Number of wavelength/velocity bins. :type wave_bins: int :param delay_bins: Number of delay time bins. :type delay_bins: int :param template: Other TF to copy all filter settings from. Will match delay, wave and velocity bins exactly. :type template: TransferFunction :param template_different_line: Is this TF going to share delay & velocity bins but have different wavelength bins? :type template_different_line: bool :param template_different_spectrum: Is this TF going to share all specified bins but be taken on photons from a different observer. :type template_different_spectrum: bool .. todo:: Consider making it impossible to apply filters after calling run(). .. py:method:: __getstate__() -> dict Removes invalid data before saving to disk. :returns: Updated internal dict, with references to external, session-specific database things, removed. :rtype: dict .. py:method:: __setstate__(state: dict) Restores the data from disk, and sets a flag to show this is a frozen TF. :param state: The unpickled object dict.. :type state: dict .. py:attribute:: _database .. py:attribute:: _session .. py:attribute:: _query .. py:attribute:: _delay_dynamic_range :value: None .. py:attribute:: _velocity :value: None .. py:attribute:: _line_list :value: None .. py:attribute:: _line_wave :value: None .. py:attribute:: _line_num :value: None .. py:attribute:: _delay_range :value: None .. py:attribute:: _continuum .. py:attribute:: _filename .. py:attribute:: _bins_wave_count :value: None .. py:attribute:: _bins_delay_count :value: None .. py:attribute:: _bins_vel :value: None .. py:attribute:: _bins_wave :value: None .. py:attribute:: _bins_delay :value: None .. py:attribute:: _emissivity :value: None .. py:attribute:: _response :value: None .. py:attribute:: _count :value: None .. py:attribute:: _wave_range :value: None .. py:attribute:: _spectrum :value: None .. py:attribute:: _unpickled :value: False .. py:method:: spectrum(number: int) -> TransferFunction Constrain the TF to photons from a specific observer :param number: Observer number from Python run :type number: int :returns: Self, so filters can be stacked :rtype: TransferFunction .. py:method:: 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 :param number: Python line number. Will vary based on data file! :type number: int :param wavelength: Wavelength of the line in angstroms :type wavelength: float :returns: Self, so filters can be stacked :rtype: TransferFunction .. py:method:: velocities(velocity: float) -> TransferFunction Constrain the TF to only photons with a range of Doppler shifts :param velocity: Maximum doppler shift velocity in m/s. Applies to both positive and negative Doppler shift :type velocity: float :returns: Self, so filters can be stacked :rtype: TransferFunction .. py:method:: wavelengths(wave_min: float, wave_max: float) -> TransferFunction Constrain the TF to only photons with a range of wavelengths :param wave_min: Minimum wavelength in angstroms :type wave_min: float :param wave_max: Maximum wavelength in angstroms :type wave_max: float :returns: Self, so filters can be stacked :rtype: TransferFunction .. py:method:: 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 :param wave_range: Array of bins to use :type wave_range: np.ndarray :returns: Self, so filters can be stacked :rtype: TransferFunction .. py:method:: 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! :param line_list: List of lines :type line_list: List[int] :returns: Self, so filters can be stacked :rtype: TransferFunction .. py:method:: delays(delay_min: float, delay_max: float, days: bool = True) -> TransferFunction The delay range that should be considered when producing the TF. :param delay_min: Minimum delay time (in seconds or days) :type delay_min: float :param delay_max: Maximum delay time (in seconds or days) :type delay_max: float :param days: Whether or not the delay range has been provided in days :type days: bool :returns: Self, so filters can be stacked :rtype: TransferFunction .. py:method:: 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 :param delay_dynamic_range: The dynamic range to be used when :type delay_dynamic_range: float :returns: Self, so filters can be stacked :rtype: TransferFunction .. py:method:: cont_scatters(scat_min: int, scat_max: Optional[int] = None) -> TransferFunction Constrain the TF to only photons that have scattered min-max times via a continuum scattering process (e.g. electron scattering). :param scat_min: Minimum number of continuum scatters :type scat_min: int :param scat_max: Maximum number of continuum scatters, if desired :type scat_max: Optional[int] :returns: Self, so filters can be stacked :rtype: TransferFunction .. py:method:: res_scatters(scat_min: int, scat_max: Optional[int] = None) -> TransferFunction Constrain the TF to only photons that have scattered min-max times via a resonant scattering process (e.g. line scattering). :param scat_min: Minimum number of resonant scatters :type scat_min: int :param scat_max: Maximum number of resonant scatters, if desired :type scat_max: Optional[int] :returns: Self, so filters can be stacked :rtype: TransferFunction .. py:method:: filter(*args) -> TransferFunction Apply a SQLalchemy filter directly to the content. :param args: The list of filter arguments :returns: Self, so filters can be stacked :rtype: TransferFunction .. py:method:: 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. :param low_state: A full, processed transfer function for a lower-luminosity system. :type low_state: TransferFunction :param high_state: A full, processed transfer function for a higher-luminosity system. :type high_state: TransferFunction :param cf_low: Correction factor for low state. Multiplier to the whole transfer function. :type cf_low: float :param cf_high: Correction factor for high state. Multiplier to the whole transfer function. :type cf_high: float :returns: Self, so plotting can be chained on. :rtype: TransferFunction .. py:method:: 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? :param response: Whether to calculate the FWHM of the transfer or response function :type response: bool :param velocity: Whether to return the FWHM in wavelength or velocity-space :type velocity: bool :returns: Full width at half maximum for the function. If the function is a doublet, this will not work properly. :rtype: float .. todo:: Catch doublets. .. py:method:: delay(response: bool = False, threshold: float = 0, bounds: float = None, days: bool = False) -> Union[float, Tuple[float, float, float]] Calculates the centroid delay for the current data :param response: Whether or not to calculate the delay from the response :type response: bool :param threshold: Exclude all bins with value < the threshold fraction of the peak value. Standard value used in the reverb papers was 0.8. :type threshold: float :param bounds: Return the fractional bounds (i.e. `bounds=0.25`, the function will return `[0.5, 0.25, 0.75]`). Not implemented. :type bounds: float :param days: Whether to return the delay in days or seconds :type days: bool :returns: Centroid delay, and lower and upper fractional bounds if bounds keyword provided :rtype: Union[float, Tuple[float, float, float]] .. todo:: Implement fractional bounds. Should just be able to call the centroid_delay function! .. py:method:: 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. :param response: Whether or not to calculate the peak transfer or response function. :type response: bool :param days: Whether to return the value in seconds or days. :type days: bool :returns: The peak delay. :rtype: float .. py:method:: 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. :param scaling_factor: 1/Number of cycles in the spectra file :type scaling_factor: float :param limit: Number of photons to limit the TF to, for testing. Recommend testing filters on a small number of photons to begin with. :type limit: int :param verbose: Whether to output exactly what the query is. :type verbose: bool :returns: Self, for chaining commands :rtype: TransferFunction .. py:method:: _return_array(array: numpy.ndarray, delay: Optional[float] = None, wave: Optional[float] = None, delay_index: Optional[int] = None) -> Union[int, float, numpy.ndarray] Internal function used by response(), emissivity() and count() :param array: Array to return value from :type array: np.ndarray :param delay: Delay to return value for. Must provide this or delay_index. :type delay: Optional[float] :param delay_index: Delay index to return value for. Must provide this or delay. :type delay_index: Optional[int] :param wave: Wavelength to return value for :type wave: Optional[float] :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. :rtype: Union[np.ndarray, float] .. todo:: Allow for only wavelength to be provided? .. py:method:: response_total() -> float Returns the total response. :returns: Total response. :rtype: float .. py:method:: delay_bins() -> numpy.ndarray Returns the range of delays covered by this TF. :returns: Array of the bin boundaries. :rtype: np.ndarray .. py:method:: response(delay: Optional[float] = None, wave: Optional[float] = None, delay_index: Optional[int] = None) -> Union[float, numpy.ndarray] Returns the responsivity in either one specific wavelength/delay bin, or all wavelength bins for a given delay. :param delay: Delay to return value for. Must provide this or delay_index. :type delay: Optional[float] :param delay_index: Delay index to return value for. Must provide this or delay. :type delay_index: Optional[int] :param wave: Wavelength to return value for. :type wave: Optional[float] :returns: Either the responsivity in one specific bin, or if wave is not specified the counts in each wavelength bin at this delay :rtype: Union[int, np.ndarray] .. todo:: Allow for only wavelength to be provided? .. py:method:: emissivity(delay: Optional[float] = None, wave: Optional[float] = None, delay_index: Optional[int] = None) -> Union[float, numpy.ndarray] Returns the emissivity in either one specific wavelength/delay bin, or all wavelength bins for a given delay. :param delay: Delay to return value for. Must provide this or delay_index. :type delay: Optional[float] :param delay_index: Delay index to return value for. Must provide this or delay. :type delay_index: Optional[int] :param wave: Wavelength to return value for. :type wave: Optional[float] :returns: Either the emissivity in one specific bin, or if wave is not specified the counts in each wavelengthin bin at this delay :rtype: Union[int, np.ndarray] .. todo:: Allow for only wavelength to be provided? .. py:method:: count(delay: Optional[float] = None, wave: Optional[float] = None, delay_index: Optional[int] = None) -> Union[int, numpy.ndarray] Returns the photon count in either one specific wavelength/delay bin, or all wavelength bins for a given delay. :param delay: Delay to return value for. Must provide this or delay_index. :type delay: Optional[float] :param delay_index: Delay index to return value for. Must provide this or delay. :type delay_index: Optional[int] :param wave: Wavelength to return value for :type wave: Optional[float] :returns: Either the count in one specific bin, or if wave is not specified the counts in each wavelength bin at this delay :rtype: Union[int, np.ndarray] .. todo:: Allow for only wavelength to be provided? .. py:method:: 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. :param response: Whether or not to return the response function data :type response: bool :param days: Whether the bin midpoints should be in seconds or days :type days: bool :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. :rtype: np.ndarray .. py:method:: 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: Optional[float] = None, format: str = '.eps', return_figure: bool = False) -> Union[TransferFunction, matplotlib.figure.Figure] Takes the data gathered by calling 'run' and outputs a plot :param log: Whether the plot should be linear or logarithmic. :type log: bool :param normalised: Whether or not to rescale the plot such that the total emissivity = 1. :type normalised: bool :param rescaled: Whether or not to rescale the plot such that the maximum emissivity = 1. :type rescaled: bool :param velocity: Whether the plot X-axis should be velocity (true) or wavelength (false). :type velocity: bool :param name: 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'. :type name: Optional[str] :param days: Whether the plot Y-axis should be in days (true) or seconds (false). :type days: bool :param response_map: Whether to plot the transfer function map or the response function. :type response_map: bool :param keplerian: 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). :type keplerian: Optional[dict] :param dynamic_range: 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. :type dynamic_range: Optional[int] :param max_delay: The optional maximum delay to plot out to. :type max_delay: Optional[float] :param rms: Whether or not the line profile panel should show the root mean squared line profile. :type rms: bool :param show: Whether or not to display the plot to screen. :type show: bool :param format: The output file format. .eps by default. :type format: str :param return_figure: If true, return the figure instead of platting it. :returns: Self, for chaining outputs :rtype: TransferFunction .. py:function:: 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! :param file_root: Root of the filename (no '.db' or '.delay_dump') :type file_root: string :param user: Username. Here in case I change to PostgreSQL :type user: string :param password: Password. Here in case I change to PostgreSQL :type password: string :param batch_size: Number of photons to stage before committing. If too low, file creation is slow. If too high, get out-of-memory errors. :type batch_size: int :returns: Connection to the database opened .. py:data:: Base Base class declared dynamically to bind to SQLalchemy .. py:class:: Spectrum Bases: :py:obj:`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. .. py:attribute:: __tablename__ :value: 'Spectra' .. py:attribute:: id .. py:attribute:: angle .. py:class:: Origin Bases: :py:obj:`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 .. py:attribute:: __tablename__ :value: 'Origin' .. py:attribute:: id .. py:attribute:: name .. py:class:: Photon Bases: :py:obj:`Base` SQLalchemy class for a photon. Why are all the properties capitalised? Changing them to lowercase as would make sense breaks backwards compatibility. # ToDo: Change to lower case. .. py:attribute:: __tablename__ :value: 'Photons' .. py:attribute:: id .. py:attribute:: Wavelength .. py:attribute:: Weight .. py:attribute:: X .. py:attribute:: Y .. py:attribute:: Z .. py:attribute:: ContinuumScatters .. py:attribute:: ResonantScatters .. py:attribute:: Delay .. py:attribute:: Spectrum .. py:attribute:: Origin .. py:attribute:: Resonance .. py:attribute:: LineResonance .. py:attribute:: Origin_matom .. py:data:: kep_sey The default Keplerian outline settings for the Seyfert model .. py:data:: kep_qso The default Keplerian outline settings for the QSO model