debrispy.asd

Functions

adapt(a_lo, a_hi, level, x1, w1, x2, w2, ...)

Perform Gauss-Legendre quadrature at two levels (lower and higher order) and adaptively refine the result until the desired tolerance is achieved or the maximum recursion depth is reached.

chunked(iterable, batch_size)

Return successive chunks of a given size from an iterable.

gl_quad(a_lo, a_hi, nodes, weights, kernel, ...)

Perform Gauss-Legendre quadrature for ASD over [a_lo, a_hi].

Classes

ASD(kernel, sigma_a)

Computes the azimuthally averaged surface density profile, ASD, from a given semi-major axis surface density Sigma(a) and kernel Phi(r, a).
class debrispy.asd.ASD(kernel, sigma_a)[source]

Bases: object

Computes the azimuthally averaged surface density profile, ASD, from a given semi-major axis surface density Sigma(a) and kernel Phi(r, a).

This class supports all kernel types.

For most use simple cases, the compute_quadvec method is the most convenient, and does not require any additional tuning of parameters, or knowledge of the integration limits.

For more complex cases, the compute_gl method is more flexible, this uses adaptive Gauss-Legendre quadrature, the user can specify the number of integration points, and the relative and absolute tolerances. Adaptive limits are supported and recommended for cases with very sharp features in the eccentricity profile.

Adaptive limits are not supported for calculations involving eccentricity distributions (only for unique eccentricity profiles).

Parameters:

  • kernel (Kernel) – A kernel object used to compute Phi(r, a).

  • sigma_a (SigmaA) – A surface density profile object defining Sigma(a).

compute_gl(r_vals=None, n_points: int = 64, tol_rel: float = 1e-08, tol_abs: float = 1e-08, max_level: int = 15, n_jobs: int = -1, show_progress: bool = True, pad: float = 0.05, rf: float = 5.0, adaptive_limits: bool = False, batch_size: int = 10, verbose: bool = True)[source]

Compute the azimuthally averaged surface density profile, ASD, using adaptive Gauss-Legendre quadrature.

Integration is performed using recursive adaptive Gauss-Legendre quadrature with optional adaptive integration bounds for efficiency (recommended for unique eccentricity profiles).

Parameters:

  • r_vals (array-like or None) – Array of radius values at which to compute the ASD. If None, defaults to 500 evenly spaced points between kernel.r_min and kernel.r_max.

  • n_points (int, optional) – Number of Gauss-Legendre points per subinterval (default: 64).

  • tol_rel (float, optional) – Relative tolerance for adaptive integration (default: 1e-8).

  • tol_abs (float, optional) – Absolute tolerance for adaptive integration (default: 1e-8).

  • max_level (int, optional) – Maximum recursion depth for adaptive integration (default: 15).

  • n_jobs (int, optional) – Number of parallel jobs to use (-1 for all available CPUs).

  • show_progress (bool, optional) – Whether to show a progress bar using tqdm (default: True).

  • pad (float, optional) – Padding added to adaptive integration bounds as a fraction of a (default: 0.05).

  • rf (float, optional) – Rescue factor: multiplier for n_points in rescue passes for suspicious results (default: 5.0).

  • adaptive_limits (bool, optional) – Whether to use adaptive integration bounds based on the eccentricity profile (default: False). Only supported for UniqueEccentricity kernels.

  • batch_size (int, optional) – Batch size for parallel integration jobs (default: 10).

  • verbose (bool, optional) – Whether to print verbose output (default: True).

Raises:

  • ValueError – If adaptive limits are enabled but the kernel does not use a unique eccentricity profile.

  • Prints –

  • ------ –

  • - Progress bar (optional) and integration diagnostics. –

  • - Warnings about any suspicious or unresolved zero results. –

Returns:

Results are stored in self._r_vals and self._sigma_r_vals.

Return type:

None

compute_quadvec(r_vals: float | ndarray[tuple[int, ...], dtype[float64]], tol_rel: float = 1e-08, tol_abs: float = 1e-08) → None[source]

Compute and cache the azimuthally averaged surface density using vectorised adaptive quadrature via scipy.integrate.quad_vec.

This method uses non-adaptive integration limits.

Parameters:

  • r_vals (float or array-like) – Radius or array of radii at which to compute ASD. The result is cached internally and accessible via self._sigma_r_vals.

  • tol_rel (float, optional) – Relative tolerance for the integration.

  • tol_abs (float, optional) – Absolute tolerance for the integration.

convolve(width: float, M: int = 2048) → ndarray[source]

Apply a Gaussian convolution to the computed ASD profile.

Parameters:

  • width (float) – The Gaussian kernel width (in the same units as r) to smooth ASD.

  • M (int, optional) – Number of uniform samples used if r-grid is non-uniform. Default is 2048.

Returns:

The convolved ASD values at the original r grid.

Return type:

np.ndarray

get_values()[source]

Return the cached ASD values.

Returns:

The cached azimuthally averaged surface density values.

Return type:

ndarray

Raises:

RuntimeError if compute_sigma_r() has not been called yet. –

integrand(a: float | ndarray[tuple[int, ...], dtype[float64]], r_vals: float | ndarray[tuple[int, ...], dtype[float64]]) → ndarray[tuple[int, ...], dtype[float64]][source]

Compute the 2D integrand over a grid of a and r values.

Parameters:

  • a (float or array) – Semi-major axis value(s) at which to evaluate the integrand.

  • r_vals (float or array) – Radial location(s) at which the azimuthally averaged profile is to be computed.

Returns:

The evaluated integrand values over all combinations of a and r_vals.

Return type:

(A, R) array

plot(a_vals: ndarray, plot_a: bool = True, x_lim: Tuple[float, float] | None = None, y_lim: Tuple[float, float] | None = None, save: bool = False, filename: str | None = None, title: str | None = None, grid_hist: bool = False) → None[source]

Plot Sigma(a) and ASD on the same figure, optionally including a histogram of the r-grid point distribution.

Parameters:

  • a_vals (array-like) – Array of semi-major axis values at which to evaluate Σ(a).

  • plot_a (bool, default=True) – Whether to plot Sigma(a) alongside ASD.

  • x_lim (tuple, optional) – Limits for the x-axis.

  • y_lim (tuple, optional) – Limits for the y-axis.

  • save (bool, default=False) – If True, saves the figure to file instead of displaying it.

  • filename (str, optional) – Name of the output file (without extension) if saving the figure.

  • title (str, optional) – Title for the plot.

  • grid_hist (bool, default=False) – If True, adds a histogram of r grid points above the main plot.

plot_convolution(width: float | None = None, M: int | None = None, overlay_sigma_r: bool = True, overlay_sigma_a: bool = False, x_lim: Tuple[float, float] | None = None, y_lim: Tuple[float, float] | None = None, save: bool = False, filename: str | None = None, title: str | None = None) → None[source]

Plot the original and convolved ASD profiles.

Parameters:

  • width (float, optional) – The Gaussian kernel width (in the same units as r) to smooth ASD. Required if convolve() has not been called yet.

  • M (int, optional) – Number of uniform samples used if r-grid is non-uniform. Required if convolve() has not been called yet.

  • overlay_sigma_r (bool, optional) – If True, overlay the original ASD profile on the convolved profile.

  • overlay_sigma_a (bool, optional) – If True, overlay the Sigma_a profile on the convolved profile.

Returns:

Displays the plot.

Return type:

None

plot_integrand(a_vals: ndarray[tuple[int, ...], dtype[float64]] | None = None, r_vals: ndarray[tuple[int, ...], dtype[float64]] | None = None, *, save: bool = False, filename: str | None = None, vmin: float | None = None, vmax: float | None = None, log: bool = True, cmap: str = 'viridis', colorbar: bool = True, r_lim: Tuple[float, float] | None = None, a_lim: Tuple[float, float] | None = None, **imshow_kwargs) → None[source]

Plot the integrand as a 2D colourmap over (r, a). Useful for diagnosing sharp features or integration difficulties.

Parameters:

  • a_vals (array-like, optional) – Semi-major axis values to evaluate the integrand over. If None, uses a default linspace.

  • r_vals (array-like, optional) – Radius values to evaluate the integrand over. If None, uses a default linspace.

  • save (bool, optional) – If True, the plot is saved to disk. Default is False.

  • filename (str, optional) – File path to save the plot. Required if save is True.

  • vmin (float, optional) – Minimum value for the colour scale. If None, uses the minimum non-zero value.

  • vmax (float, optional) – Maximum value for the colour scale. If None, uses the maximum value.

  • log (bool, optional) – If True, applies logarithmic colour scaling (ignoring non-positive values). Default is True.

  • cmap (str, optional) – Colormap to use for the image. Default is “viridis”.

  • colorbar (bool, optional) – If True, adds a colourbar to the plot. Default is True.

  • r_lim (tuple of float, optional) – Limits for the x-axis (r-axis).

  • a_lim (tuple of float, optional) – Limits for the y-axis (a-axis).

  • **imshow_kwargs (dict) – Additional keyword arguments passed to imshow().

Raises:

ValueError – If save=True but filename is not provided.

refine(curvature_factor: float = 1.0, max_rounds: int = 3, subdiv: int = 2, n_jobs: int = -1, show_progress: bool = False, n_points: int = 64, tol_rel: float = 1e-08, tol_abs: float = 1e-08, max_level: int = 25, adaptive_limits: bool = False, rf: float = 10.0, pad: float = 0.05, batch_size: int = 10)[source]

Refine the existing ASD grid by adding intermediate points based on curvature.

This method uses a curvature-driven refinement strategy: it identifies regions where the ASD changes rapidly and adds more sampling points in those regions. Optionally, it applies adaptive integration bounds based on the eccentricity profile and performs a rescue pass on suspicious zeros.

This uses the same adaptive Gauss-Legendre quadrature as the compute_sigma_r method. Currently, this is only available for unique eccentricity profiles.

Parameters:

  • curvature_factor (float, optional) – Threshold multiplier for triggering refinement based on curvature.

  • max_rounds (int, optional) – Maximum number of refinement rounds to perform. Each round may add new points.

  • subdiv (int, optional) – Number of subdivisions to insert between points with high curvature. For example, subdiv=2 inserts one point between each flagged pair.

  • n_jobs (int, optional) – Number of parallel jobs to use for integration (-1 = all available CPUs).

  • show_progress (bool, optional) – Whether to show a progress bar during integration.

  • n_points (int, optional) – Number of Gauss-Legendre nodes to use in initial integration.

  • tol_rel (float, optional) – Relative tolerance for the adaptive integrator.

  • tol_abs (float, optional) – Absolute tolerance for the adaptive integrator.

  • max_level (int, optional) – Maximum recursion depth for the adaptive Gauss-Legendre integrator.

  • adaptive_limits (bool, optional) – If True, use adaptive (localised) integration bounds based on the eccentricity profile. Requires that the kernel’s eccentricity profile is UniqueEccentricity.

  • rf (float, optional) – Rescue factor: multiply n_points by this value during the rescue pass to resolve suspicious zeros.

  • pad (float, optional) – Padding to apply to adaptive integration bounds (in units of semi-major axis).

  • batch_size (int, optional) – Batch size for parallel job distribution.

Raises:

ValueError – If adaptive_limits=True but the kernel does not have a UniqueEccentricity profile.

Returns:

Updates the internal self._r_vals and self._sigma_r_vals arrays with the refined grid.

Return type:

None

debrispy.asd.adapt(a_lo: float, a_hi: float, level: int, x1: ndarray[tuple[int, ...], dtype[float64]], w1: ndarray[tuple[int, ...], dtype[float64]], x2: ndarray[tuple[int, ...], dtype[float64]], w2: ndarray[tuple[int, ...], dtype[float64]], tol_rel: float, tol_abs: float, max_level: int, kernel, sigma_a, r_j: float) → float[source]

Perform Gauss-Legendre quadrature at two levels (lower and higher order) and adaptively refine the result until the desired tolerance is achieved or the maximum recursion depth is reached.

Parameters:

  • a_lo (float) – Lower limit of integration.

  • a_hi (float) – Upper limit of integration.

  • level (int) – Current recursion depth.

  • x1 ((N,) arrays) – Lower-order Gauss-Legendre nodes and weights.

  • w1 ((N,) arrays) – Lower-order Gauss-Legendre nodes and weights.

  • x2 ((N,) arrays) – Higher-order Gauss-Legendre nodes and weights.

  • w2 ((N,) arrays) – Higher-order Gauss-Legendre nodes and weights.

  • tol_rel (float) – Relative tolerance.

  • tol_abs (float) – Absolute tolerance.

  • max_level (int) – Maximum recursion depth.

  • kernel (Kernel) – Kernel object.

  • sigma_a (SigmaA) – Surface density profile.

  • r_j (float) – Radius value at which the integral is being computed.

Returns:

Integral value of ASD.

Return type:

float

debrispy.asd.chunked(iterable: Iterable[T], batch_size: int) → Iterator[List[T]][source]

Return successive chunks of a given size from an iterable.

Parameters:

  • iterable (Iterable[T]) – Input iterable to be split into chunks.

  • batch_size (int) – Maximum size of each chunk.

Returns:

Batches of the original iterable as lists of size at most batch_size.

Return type:

Iterator[List[T]]

debrispy.asd.gl_quad(a_lo: float, a_hi: float, nodes: ndarray[tuple[int, ...], dtype[float64]], weights: ndarray[tuple[int, ...], dtype[float64]], kernel, sigma_a, r_j: float) → float[source]

Perform Gauss-Legendre quadrature for ASD over [a_lo, a_hi].

Parameters:

  • a_lo (float) – Lower integration bound.

  • a_hi (float) – Upper integration bound.

  • nodes ((N,) array) – Gauss-Legendre nodes in [-1, 1].

  • weights ((N,) array) – Corresponding Gauss-Legendre weights.

  • kernel (Kernel object) – Phi(r, a) evaluator.

  • sigma_a (SigmaA object) – Sigma(a) surface density profile.

  • r_j (float) – Evaluation point in radius.

Returns:

Integral value of ASD.

Return type:

float