dased.criteria¶

Optimization criteria for DAS layout design.

This package exposes classes such as EIGCriterion and RaySensitivity used by the optimization routines in optimisation.

Classes¶

`EIGCriterion`	Expected Information Gain (EIG) criterion for Bayesian experimental design of DAS layouts.
`RaySensitivity`	Ray-based sensitivity criterion for optimal DAS layout design.

Functions¶

get_gaussian_prior

Create a Gaussian prior covariance matrix for a 2D grid using vectorized operations.

Descriptions¶

class dased.criteria.EIGCriterion(samples, data_likelihood, eig_method='NMC', random_seed=0, downsample=None, **kwargs)¶

Expected Information Gain (EIG) criterion for Bayesian experimental design of DAS layouts.

This class evaluates DAS layout designs based on the Expected Information Gain (EIG), which quantifies how much information a design is expected to provide about model parameters given prior knowledge and a data likelihood function. The implementation relies on the geobed package for EIG computation.

Parameters:

samples (ndarray | Tensor) –
Prior samples of model parameters as a 2D array/tensor with shape (N_samples, N_parameters). Each row represents one sample from the prior distribution. Can be a ndarray

or torch.Tensor. The implementation converts inputs to torch.Tensor internally for compatibility with geobed.
data_likelihood (callable) –
Callable function that computes the likelihood of observed data given model parameters and experimental design. Should have signature:
data_likelihood(model_params, design_gdf, data) -> distribution

Where: - model_params: Current parameter values from prior samples - design_gdf: geopandas.GeoDataFrame representing the

DAS layout (see get_gdf)
- data: Observed or simulated data
- distribution: A torch.distributions.Distribution
  instance representing the likelihood of the data given the parameters and design.
eig_method (str) –
Method for EIG calculation. Supported methods depend on geobed implementation. Common options include:
- "NMC": Nested Monte Carlo (default)
- "DN": \(\mathrm{D}_N\)-Method
Defaults to “NMC”.
random_seed (int) – Random seed for reproducible EIG calculations. Ensures consistent results across multiple runs with the same inputs. Defaults to 0.
downsample (int) –
Optional downsampling factor for reducing computational cost. If specified, reduces the number of channels in the layout by selecting representative channels from blocks of size downsample. Useful when full resolution is not required.

If None, uses all channels. If > 1, applies downsampling. Defaults to None.
**kwargs (Any) –
Additional keyword arguments passed to the EIG calculation method. Common options include:
- N: Number of samples to use (default: min(1000, n_samples))
- reuse_M: Whether to reuse samples in the inner loop of the NMC method for efficiency (NMC only)
See geobed documentation for method-specific parameters.

Raises:

ImportError – If geobed package is not installed.
TypeError – If samples is not a numpy array or torch tensor, or if data_likelihood is not callable.
ValueError – If samples is not 2D, has fewer than required samples for the chosen method, or if EIG method parameters are invalid.

Examples

Basic EIG criterion setup:

>>> import numpy as np
>>> from geopandas import GeoDataFrame
>>>
>>> # Generate prior samples (e.g., velocity model parameters) :class:`~numpy.ndarray`
>>> prior_samples = np.random.normal(3000, 300, (1000, 2))  # :class:`~numpy.ndarray`
>>>
>>> # Define likelihood function
>>> def likelihood_func(params, design_gdf, data):
...     # Example likelihood: Gaussian with mean based on params and design
...     mean = f(params, design_gdf)  # User-defined function
...     return torch.distributions.Normal(loc=mean, scale=1.0)  # :class:`torch.distributions.Normal`
>>>
>>> # Create EIG criterion
>>> eig_criterion = EIGCriterion(
...     samples=prior_samples,
...     data_likelihood=likelihood_func,
...     eig_method="NMC"
... )

Evaluating a DAS layout:

>>> from dased.layout import DASLayout
>>>
>>> # Create layout (knots as a NumPy array :class:`~numpy.ndarray`)
>>> knots = np.array([[0, 0], [1000, 0]])  # :class:`~numpy.ndarray`
>>> layout = DASLayout(knots, spacing=10.0)
>>>
>>> # Calculate EIG
>>> eig_value = eig_criterion(layout)
>>> print(f"Expected Information Gain: {eig_value:.3f}")

__call__(design)¶

Calculate the Expected Information Gain for a given DAS layout design.

This method evaluates how much information the specified layout is expected to provide about the model parameters, considering the prior uncertainty and the data likelihood function.

Parameters:: design (DASLayout | GeoDataFrame) – DAS layout design to evaluate. Can be either:

:param - DASLayout: Layout object containing channel positions and properties :param - geopandas.GeoDataFrame: Direct GeoDataFrame representation of the layout

The design must contain valid geometry (Point objects for channels) and any required properties for the likelihood function.

Returns:

Expected Information Gain value as a float. Higher values indicate designs that are expected to provide more information about the model parameters.

Raises:

TypeError – If design is not a DASLayout or GeoDataFrame
ValueError – If design lacks required geometry or has invalid structure, or if downsampling parameters result in no remaining channels

Return type:

float

Examples

Evaluate a simple linear layout:

>>> import numpy as np
>>> from dased.layout import DASLayout
>>>
>>> # Create layout (knots as a NumPy array :class:`~numpy.ndarray`)
>>> knots = np.array([[0, 0], [1000, 0]])  # :class:`~numpy.ndarray`
>>> layout = DASLayout(knots, spacing=25.0)
>>>
>>> # Calculate EIG (assuming criterion is already configured)
>>> eig_value = eig_criterion(layout)
>>> print(f"EIG: {eig_value:.3f}")

Evaluate with downsampling:

>>> # Create criterion with downsampling
>>> eig_criterion_ds = EIGCriterion(
...     samples=prior_samples,
...     data_likelihood=likelihood_func,
...     downsample=4  # Use every 4th channel
... )
>>> eig_value_ds = eig_criterion_ds(layout)

class dased.criteria.RaySensitivity(x_range, y_range, data_type, n_points=70, reference_distance=1000.0, snr_threshold=1.5, sources=None, roi=None, criterion='D', criterion_kwargs=None)¶

Ray-based sensitivity criterion for optimal DAS layout design.

This class evaluates DAS layouts using straight-ray path assumptions to compute sensitivity-based optimality criteria. The method constructs a linearized sensitivity matrix relating model parameters (for example, seismic velocity) to data observations, then computes A-optimality, D-optimality, or RER measures.

Type hints and documentation use Sphinx roles for common external types such as ndarray, pandas.DataFrame, and geopandas.GeoDataFrame. Internal references (for example DASLayout) use fully qualified references.

The criterion supports both source-to-receiver and inter-channel (ambient noise) ray configurations, with Signal-to-Noise Ratio (SNR) filtering based on:

Geometric spreading (distance decay)
Wave type directional sensitivity (Rayleigh vs Love waves)
Channel-specific signal strength and coupling factors
User-defined SNR thresholds

Parameters:

x_range (Tuple[float, float]) – Spatial extent in x-direction as (min, max) tuple in meters. Defines the horizontal bounds of the model parametrization grid.
y_range (Tuple[float, float]) – Spatial extent in y-direction as (min, max) tuple in meters. Defines the vertical bounds of the model parametrization grid.
data_type (str) –
Wave type for directional sensitivity calculations. Supported types:
- "rayleigh": Rayleigh waves with \(\cos^2(\theta)\) directional scaling
- "love": Love waves with \(|\sin(2\theta)|\) directional scaling
- "3C": Three-component (no directional scaling)
Where \(\theta\) is the angle between ray direction and channel orientation.
n_points (int | Tuple[int, int]) –
Grid resolution for model parametrization. Can be:
- int: Creates square grid with n_points \(\times\) n_points cells
- tuple: Specifies (nx, ny) for rectangular grids
Higher resolution increases computational cost but improves spatial detail. Defaults to 70.
reference_distance (float) – Reference distance in meters for SNR normalization. Distance at which SNR = 1.0 for signal strength calculation. Used in geometric spreading: SNR \(\propto\) (distance/reference_distance):math:^{-1}. Defaults to 1000.0 meters.
snr_threshold (float) – Minimum SNR threshold for ray inclusion. Defaults to 1.5.
sources (ndarray | None) –
Optional[ndarray] External source locations for source-to-receiver ray paths. If provided, should be a 2D array with shape (M, >=2) containing (x, y) coordinates, optionally (x, y, z). If None, uses inter-channel rays for ambient noise analysis.

Format: [[x1, y1], [x2, y2], ..., [xM, yM]]. Units should match x_range and y_range (typically meters). Defaults to None.
roi (Any | None) –
Region of Interest for constraining the model space. Options:
- None: Use entire grid (default)
- list/tuple: Bounding box as [xmin, xmax, ymin, ymax]
- shapely.geometry.Polygon/shapely.geometry.MultiPolygon: Shapely geometry objects
Only grid cells within the ROI contribute to the optimization criterion. Useful for focusing on area(s) of interest. Defaults to None.
criterion (str) –
Optimality criterion for layout evaluation. Supported options:
- "A": A-optimality (trace of covariance matrix)
- "D": D-optimality (determinant of covariance matrix)
- "RER": RER Criterion
A-optimality minimizes average parameter uncertainty, D-optimality minimizes uncertainty volume, RER optimises number of eigenvalues above threshold. Defaults to “D”.
criterion_kwargs (dict | None) – Additional parameters for criterion calculation. Dictionary passed to the optimality computation method. Specific options depend on the chosen criterion. Defaults to None (empty dictionary).

Raises:

ImportError – If ttcrpy package is not installed (required for ray tracing).
ValueError – If data_type or criterion are not in supported lists, or if sources array has invalid dimensions.

Examples

Basic D-optimality criterion for Rayleigh waves:

>>> criterion = RaySensitivity(
...     x_range=(0, 5000),     # 5 km extent
...     y_range=(0, 3000),     # 3 km extent
...     data_type="rayleigh",
...     n_points=80,           # 80x80 grid
...     snr_threshold=1.5
... )

Source-receiver configuration with A-optimality:

>>> import numpy as np
>>> sources = np.array([[100, 100], [900, 100], [500, 500]])  # :class:`~numpy.ndarray`
>>> criterion = RaySensitivity(
...     x_range=(0, 1000),
...     y_range=(0, 600),
...     data_type="love",
...     sources=sources,
...     criterion="A"
... )

Using ROI to focus on specific area:

>>> from shapely.geometry import Polygon
>>> roi_poly = Polygon([(200, 200), (800, 200), (800, 400), (200, 400)])  # :class:`shapely.geometry.Polygon`
>>> criterion = RaySensitivity(
...     x_range=(0, 1000),
...     y_range=(0, 600),
...     data_type="rayleigh",
...     roi=roi_poly,
...     criterion="RER"
... )

__call__(design)¶

Calculate the optimality criterion value for a given DAS layout design.

This method evaluates the specified design by constructing a sensitivity matrix based on ray paths and computing the selected optimality criterion (A-optimality, D-optimality, or RER).

Parameters:

design (DASLayout | GeoDataFrame) –

DAS layout design to evaluate. Can be either:

DASLayout: Layout object with channel positions and properties
GeoDataFrame: Direct GeoDataFrame representation with Point geometries

Must contain columns ‘u_x’ and ‘u_y’ for channel orientations. Optional columns include ‘signal_strength’, ‘K’, and data-type specific ‘K_rayleigh’ or ‘K_love’ for SNR calculations.

Returns:

A-optimality: Higher values indicate better designs (minimized average uncertainty)
D-optimality: Higher values indicate better designs (maximized information content)
RER: Higher values indicate better RER Criterion

Returns 0.0 if no valid rays meet the SNR threshold.

Return type:

Optimality criterion value as a float. Interpretation depends on criterion

Raises:

TypeError – If design is not a DASLayout or GeoDataFrame
ValueError – If design lacks required geometry or orientation columns

Examples

Evaluate a linear array for ambient noise tomography:

>>> import numpy as np
>>> from dased.layout import DASLayout
>>>
>>> # Create simple linear layout
>>> knots = np.array([[0, 0], [2000, 0]])  # :class:`~numpy.ndarray`
>>> layout = DASLayout(knots, spacing=50.0)
>>>
>>> # Evaluate with A-optimality
>>> score = criterion(layout)  # assuming criterion is configured
>>> print(f"A-optimality score: {score:.3f}")

Compare different layout configurations:

>>> # L-shaped layout
>>> knots_L = np.array([[0, 0], [1000, 0], [1000, 1000]])  # :class:`~numpy.ndarray`
>>> layout_L = DASLayout(knots_L, spacing=50.0)
>>>
>>> # Circular layout
>>> theta = np.linspace(0, 2*np.pi, 21)[:-1]  # 20 points
>>> knots_circle = 500 * np.column_stack([np.cos(theta), np.sin(theta)])  # :class:`~numpy.ndarray`
>>> layout_circle = DASLayout(knots_circle, spacing=50.0)
>>>
>>> score_L = criterion(layout_L)
>>> score_circle = criterion(layout_circle)
>>> print(f"L-shape: {score_L:.3f}, Circle: {score_circle:.3f}")

plot(design, title=None, ax=None, show_stats=True, plot_sources=True, plot_roi_boundary=True, log_scale=False, **kwargs)¶

Plots the sensitivity map and related info.

Parameters:

design (DASLayout) – DASLayout object.
title (str | None) – Plot title.
ax (Axes | None) – Axes to plot on. If None, creates a new figure/axes.
show_stats (bool) – Display statistics text box.
plot_sources (bool) – Plot source markers if sources were provided.
plot_roi_boundary (bool) – Plot the ROI boundary if ROI was provided.
log_scale (bool) – If True, apply log1p transformation to sensitivity values.
**kwargs – Additional arguments passed to ax.pcolormesh.

Returns:

(fig, ax) The figure and axes objects, or (None, None) on failure.

Return type:

tuple

get_eigenvalue_spectrum(design, normalise=True)¶

Computes the normalized eigenvalue spectrum of the sensitivity matrix for the given design.

Parameters:: design (DASLayout) – DASLayout object.
Returns:: Normalized eigenvalues in descending order.
Return type:: np.ndarray

plot_eigenvalue_spectrum(design, title=None, ax=None, **kwargs)¶

Plots the eigenvalue spectrum of the sensitivity matrix.

Parameters:

design (DASLayout) – DASLayout object.
title (str | None) – Plot title.
ax (Axes | None) – Axes to plot on. If None, creates a new figure/axes.
**kwargs – Additional arguments passed to ax.plot.

Returns:

(fig, ax) The figure and axes objects, or (None, None) on failure.

Return type:

tuple

inverse(design, m_true, m_prior, sigma_d, correlation_length, regularization_weight)¶

Performs a simple linear inversion to estimate the slowness model.

Parameters:

design – DASLayout or GeoDataFrame for the layout.
m_true – True slowness model (2D array).
m_prior – Prior slowness model (2D array).
sigma_d – Data standard deviation (float).
correlation_length – Prior correlation length (float).
regularization_weight – Prior regularization weight (float).

Returns:

Estimated slowness model (2D array).

plot_checkerboard(design, background_velocity, perturbation, vmin, vmax, sigma_d, correlation_length, regularization_weight, ax=None, checkerboard_size=4, mask_outside_roi=True, show_metrics=True, **kwargs)¶

Plots the true and estimated checkerboard velocity pattern using the provided design.

Parameters:

design – DASLayout or GeoDataFrame for the layout.
ax (Axes | None) – Optional matplotlib axes to plot on. If provided, only the estimated velocity is plotted.
checkerboard_size (int) – Size of each block in the checkerboard.
background_velocity (float) – Background velocity (m/s).
perturbation (float) – Relative perturbation (e.g., 0.05 for 5%).
vmin (float) – Minimum velocity for color scale.
vmax (float) – Maximum velocity for color scale.
sigma_d (float) – Data standard deviation.
correlation_length (float) – Prior correlation length.
regularization_weight (float) – Prior regularization weight.
mask_outside_roi (bool) – If True, set values outside ROI to NaN.
show_metrics (bool) – If True, display resolution metrics (AR, CC, RMSE) on the plot.
**kwargs – Additional arguments for pcolormesh.

Returns:

The matplotlib axes with the plot.

dased.criteria.get_gaussian_prior(smoothing_length, nx, ny, grid_spacing)¶

Create a Gaussian prior covariance matrix for a 2D grid using vectorized operations.

This function generates a sparse covariance matrix that encodes Gaussian spatial correlation between grid cells. The covariance decreases exponentially with distance according to the specified smoothing length scale.

Parameters:

smoothing_length (float) – Characteristic length scale for Gaussian smoothing in meters. Controls the spatial correlation range - larger values create smoother spatial fields with longer-range correlations.
nx (int) – Number of grid cells in the x-direction (horizontal).
ny (int) – Number of grid cells in the y-direction (vertical).
grid_spacing (Tuple[float, float]) – Tuple containing grid spacing in (x, y) directions in meters. Format: (dx, dy) where dx and dy are the cell dimensions.

Returns:

Sparse CSC (Compressed Sparse Column) matrix of shape (N, N) where N = nx * ny. Matrix element (i, j) represents the covariance between grid cells i and j. Diagonal elements are normalized to 1.0, off-diagonal elements decay exponentially with spatial separation.

Return type:

csr_matrix

Examples

Create prior for a 50x50 grid with 100m spacing:

>>> grid_spacing = (100.0, 100.0)  # 100m x 100m cells
>>> smoothing_length = 500.0       # 500m correlation length
>>> prior_cov = get_gaussian_prior(smoothing_length, 50, 50, grid_spacing)
>>> print(f"Prior matrix shape: {prior_cov.shape}")
>>> print(f"Matrix density: {prior_cov.nnz / (50*50)**2:.3f}")

Use with different x/y spacing:

>>> # Rectangular cells: 200m x 100m
>>> grid_spacing = (200.0, 100.0)
>>> prior_cov = get_gaussian_prior(300.0, 25, 50, grid_spacing)