PDE models

sies.pde.conductivity

The conductivity problem in free space.

The potential \(u\) solves

$$ \nabla \cdot (\rho_D \nabla u) = 0 \quad \text{in } \mathbb{R}^2 \setminus {x_s}, \qquad u - G(\cdot - x_s) = O(|x|^{-1}), $$ where \(\rho_D = 1 + \sum_l (k_l - 1) \chi_{D_l}\) and \(k_l = \sigma_l + i \omega \epsilon_l\). The Multi-Static Response (MSR) matrix collects the perturbations \(u - G\) measured at the receivers; it is simulated here through a boundary integral representation, and inverted for the CGPT of the inclusions.

Reference: Ammari et al., Target identification using dictionary matching of generalized polarization tensors, FoCM (2014).

MSRData dataclass

MSRData(msr, freqs, msr_noisy=list(), noise_sigma=list())

Multi-static response data, possibly at several frequencies.

Attributes:

Name	Type	Description
msr	list of ndarray	One MSR matrix of shape (nb_sources, nb_receivers_per_group) per frequency. Entry (s, r) is the field perturbation produced by source s and measured at the r-th receiver of the source's group (for single-group configurations, simply the r-th receiver).
freqs	list of float	Working frequencies.
msr_noisy	list of ndarray	Noisy version of msr (empty until noise is added).
noise_sigma	list of float	Standard deviation of the noise added at each frequency.

ReconstructionResult dataclass

ReconstructionResult(cgpt, residual, relative_residual, source_matrix, receiver_matrix)

Result of a CGPT reconstruction.

Attributes:

Name	Type	Description
cgpt	list of ndarray	Reconstructed CGPT matrix at each frequency.
residual	list of float	Norm of the data misfit at each frequency.
relative_residual	list of float	Data misfit relative to the data norm.
source_matrix	ndarray	Linear operator As acting on the source side.
receiver_matrix	ndarray	Linear operator Ar acting on the receiver side; the forward model is MSR = As @ CGPT @ Ar.T.

ConductivityR2

ConductivityR2(inclusions, cnd, pmtt, cfg)

Conductivity problem with small inclusions in free space.

Parameters:

Name	Type	Description	Default
inclusions	C2Boundary or list of C2Boundary	Mutually disjoint inclusions, all discretized with the same number of boundary points.	required
cnd	array_like	Conductivity of each inclusion (positive, different from one).	required
pmtt	array_like	Permittivity of each inclusion (nonnegative).	required
cfg	AcquisitionConfig	Geometry of sources and receivers.	required

Attributes:

Name	Type	Description
inclusions	list of C2Boundary	The inclusions.
cnd, pmtt	ndarray	Material constants.
cfg	AcquisitionConfig	The acquisition configuration.

Source code in src/sies/pde/conductivity.py

def __init__(
    self,
    inclusions: list[C2Boundary] | C2Boundary,
    cnd: NDArray,
    pmtt: NDArray,
    cfg: AcquisitionConfig,
):
    if isinstance(inclusions, C2Boundary):
        inclusions = [inclusions]
    nb_points = {incl.nb_points for incl in inclusions}
    if len(nb_points) != 1:
        raise ValueError("All inclusions must have the same number of boundary points.")

    cnd = np.atleast_1d(np.asarray(cnd, dtype=float))
    pmtt = np.atleast_1d(np.asarray(pmtt, dtype=float))
    if len(cnd) != len(inclusions) or len(pmtt) != len(inclusions):
        raise ValueError("Conductivity and permittivity must be given for each inclusion.")
    if np.any(cnd == 1) or np.any(cnd < 0):
        raise ValueError("Conductivity must be nonnegative and different from 1.")
    if np.any(pmtt < 0):
        raise ValueError("Permittivity must be nonnegative.")

    self.inclusions = inclusions
    self.cnd = cnd
    self.pmtt = pmtt
    self.cfg = cfg

    # Frequency-independent quantities, precomputed once.
    self._block_matrix = make_block_matrix(inclusions)
    self._dgdn = self._compute_dgdn()

simulate_data

simulate_data(freqs=0.0)

Simulate the MSR matrices at the given frequencies.

The perturbation measured by a receiver is represented as a sum of single layer potentials, \(u - G = \sum_l S_{D_l}[\phi_l]\).

Parameters:

Name	Type	Description	Default
freqs	float or list of float	Working frequencies.	0.0

Returns:

Type	Description
MSRData	The simulated data.

Source code in src/sies/pde/conductivity.py

def simulate_data(self, freqs: float | list[float] = 0.0) -> MSRData:
    r"""Simulate the MSR matrices at the given frequencies.

    The perturbation measured by a receiver is represented as a sum
    of single layer potentials,
    $u - G = \sum_l S_{D_l}[\phi_l]$.

    Parameters
    ----------
    freqs : float or list of float, default 0.0
        Working frequencies.

    Returns
    -------
    MSRData
        The simulated data.
    """
    freqs = np.atleast_1d(np.asarray(freqs, dtype=float))

    msr_list = []
    for freq in freqs:
        densities = self._compute_densities(freq)
        dtype = complex if any(np.iscomplexobj(p) for p in densities) else float
        msr = np.zeros((self.cfg.nb_sources, self.cfg.nb_receivers_per_group), dtype=dtype)
        for incl, phi in zip(self.inclusions, densities, strict=True):
            for s in range(self.cfg.nb_sources):
                rcv = self.cfg.receivers_of_source(s)
                msr[s] += SingleLayer.evaluate(incl, phi[:, s], rcv)
        msr_list.append(msr)

    return MSRData(msr=msr_list, freqs=list(freqs))

add_white_noise staticmethod

add_white_noise(data, level, rng=None)

Add white noise to simulated MSR data.

The real and imaginary parts are corrupted independently, source by source.

Parameters:

Name	Type	Description	Default
data	MSRData	Simulated data.	required
level	float	Noise level (e.g. 0.1 for 10% noise).	required
rng	Generator	Random generator, for reproducibility.	None

Returns:

Type	Description
MSRData	A new data object with the msr_noisy and noise_sigma fields filled in.

Source code in src/sies/pde/conductivity.py

@staticmethod
def add_white_noise(
    data: MSRData, level: float, rng: np.random.Generator | None = None
) -> MSRData:
    """Add white noise to simulated MSR data.

    The real and imaginary parts are corrupted independently, source
    by source.

    Parameters
    ----------
    data : MSRData
        Simulated data.
    level : float
        Noise level (e.g. ``0.1`` for 10% noise).
    rng : numpy.random.Generator, optional
        Random generator, for reproducibility.

    Returns
    -------
    MSRData
        A new data object with the `msr_noisy` and `noise_sigma`
        fields filled in.
    """
    rng = rng or np.random.default_rng()
    noisy, sigmas = [], []
    for msr in data.msr:
        real, sigma_r = add_white_noise(msr.real, level, rng)
        imag, sigma_i = add_white_noise(msr.imag, level, rng)
        noisy.append(real + 1j * imag if np.iscomplexobj(msr) else real)
        sigmas.append(abs(sigma_r + 1j * sigma_i))
    return MSRData(msr=data.msr, freqs=data.freqs, msr_noisy=noisy, noise_sigma=sigmas)

make_matrix_a staticmethod

make_matrix_a(points, center, order)

Acquisition matrix of the linearized CGPT forward model.

Row i contains the coefficients \([\cos(m\theta_i), \sin(m\theta_i)] / (2\pi m R_i^m)\) for m = 1..order, where \((R_i, \theta_i)\) are the polar coordinates of the i-th point relative to center.

Parameters:

Name	Type	Description	Default
points	(ndarray, shape(2, n))	Coordinates of the sources or receivers.	required
center	(ndarray, shape(2))	Reference center.	required
order	int	Maximum CGPT order.	required

Returns:

Type	Description
(ndarray, shape(n, 2 * order))	The acquisition matrix.

Source code in src/sies/pde/conductivity.py

@staticmethod
def make_matrix_a(points: NDArray, center: NDArray, order: int) -> NDArray:
    r"""Acquisition matrix of the linearized CGPT forward model.

    Row ``i`` contains the coefficients
    $[\cos(m\theta_i), \sin(m\theta_i)] / (2\pi m R_i^m)$
    for ``m = 1..order``, where $(R_i, \theta_i)$ are the
    polar coordinates of the ``i``-th point relative to `center`.

    Parameters
    ----------
    points : ndarray, shape (2, n)
        Coordinates of the sources or receivers.
    center : ndarray, shape (2,)
        Reference center.
    order : int
        Maximum CGPT order.

    Returns
    -------
    ndarray, shape (n, 2 * order)
        The acquisition matrix.
    """
    delta = points - np.asarray(center, dtype=float).reshape(2, 1)
    radius = np.linalg.norm(delta, axis=0)
    angle = np.arctan2(delta[1], delta[0])

    m = np.arange(1, order + 1)
    weights = 1 / (2 * np.pi * m * radius[:, np.newaxis] ** m)
    matrix = np.empty((points.shape[1], 2 * order))
    matrix[:, 0::2] = np.cos(angle[:, np.newaxis] * m) * weights
    matrix[:, 1::2] = np.sin(angle[:, np.newaxis] * m) * weights
    return matrix

make_linear_operator

make_linear_operator(order)

Build the matrices of the forward model MSR = As CGPT Ar^T.

Only single-group (full or sparse view) configurations are supported.

Parameters:

Name	Type	Description	Default
order	int	Maximum CGPT order.	required

Returns:

Name	Type	Description
source_matrix	(ndarray, shape(nb_sources, 2 * order))	Matrix As acting on the source side.
receiver_matrix	(ndarray, shape(nb_receivers, 2 * order))	Matrix Ar acting on the receiver side.

Source code in src/sies/pde/conductivity.py

def make_linear_operator(self, order: int) -> tuple[NDArray, NDArray]:
    """Build the matrices of the forward model ``MSR = As CGPT Ar^T``.

    Only single-group (full or sparse view) configurations are
    supported.

    Parameters
    ----------
    order : int
        Maximum CGPT order.

    Returns
    -------
    source_matrix : ndarray, shape (nb_sources, 2 * order)
        Matrix ``As`` acting on the source side.
    receiver_matrix : ndarray, shape (nb_receivers, 2 * order)
        Matrix ``Ar`` acting on the receiver side.
    """
    if self.cfg.nb_groups != 1:
        raise NotImplementedError("Grouped (limited-view) operators are not supported.")
    src_matrix = self.make_matrix_a(self.cfg.all_sources, self.cfg.center, order)
    rcv_matrix = self.make_matrix_a(self.cfg.all_receivers, self.cfg.center, order)
    return src_matrix, rcv_matrix

reconstruct_cgpt

reconstruct_cgpt(msr, order)

Reconstruct the CGPT from MSR data by least squares.

Parameters:

Name	Type	Description	Default
msr	ndarray or list of ndarray	MSR matrix (or one matrix per frequency).	required
order	int	Maximum order of the reconstruction.	required

Returns:

Type	Description
ReconstructionResult	The reconstructed CGPT and residuals.

Source code in src/sies/pde/conductivity.py

def reconstruct_cgpt(self, msr: NDArray | list[NDArray], order: int) -> ReconstructionResult:
    """Reconstruct the CGPT from MSR data by least squares.

    Parameters
    ----------
    msr : ndarray or list of ndarray
        MSR matrix (or one matrix per frequency).
    order : int
        Maximum order of the reconstruction.

    Returns
    -------
    ReconstructionResult
        The reconstructed CGPT and residuals.
    """
    src_matrix, rcv_matrix = self.make_linear_operator(order)
    return self._invert(msr, src_matrix, rcv_matrix)

reconstruct_cgpt_analytic

reconstruct_cgpt_analytic(msr, order)

Reconstruct the CGPT with the closed-form least-squares inverse.

For equispaced full-view circular configurations the acquisition matrices satisfy Cs' Cs = Ns / 2 I, which yields the explicit solution

\[M = \frac{4}{N_s N_r} D_s^{-1} C_s^T \, \mathrm{MSR} \, C_r D_r^{-1}.\]

Parameters:

Name	Type	Description	Default
msr	ndarray or list of ndarray	MSR matrix (or one matrix per frequency).	required
order	int	Maximum order of the reconstruction. Internally capped at (min(Ns, Nr) - 1) / 2, the maximum resolvable order.	required

Returns:

Type	Description
ReconstructionResult	The reconstructed CGPT and residuals.

Source code in src/sies/pde/conductivity.py

def reconstruct_cgpt_analytic(
    self, msr: NDArray | list[NDArray], order: int
) -> ReconstructionResult:
    r"""Reconstruct the CGPT with the closed-form least-squares inverse.

    For equispaced full-view circular configurations the acquisition
    matrices satisfy ``Cs' Cs = Ns / 2 I``, which yields the explicit
    solution

    $$M = \frac{4}{N_s N_r} D_s^{-1} C_s^T \, \mathrm{MSR} \, C_r D_r^{-1}.$$

    Parameters
    ----------
    msr : ndarray or list of ndarray
        MSR matrix (or one matrix per frequency).
    order : int
        Maximum order of the reconstruction. Internally capped at
        ``(min(Ns, Nr) - 1) / 2``, the maximum resolvable order.

    Returns
    -------
    ReconstructionResult
        The reconstructed CGPT and residuals.
    """
    cfg = self.cfg
    if not (isinstance(cfg, Concentric) and cfg.equispaced):
        raise ValueError(
            "Analytic reconstruction requires an equispaced concentric configuration."
        )

    nb_src, nb_rcv = cfg.nb_sources, cfg.nb_receivers
    order = min(order, (nb_src - 1) // 2, (nb_rcv - 1) // 2)

    cs, ds = _make_matrix_cd(nb_src, cfg.radius_src, order)
    cr, dr = _make_matrix_cd(nb_rcv, cfg.radius_rcv, order)

    msr_list = [msr] if not isinstance(msr, list) else msr
    ids = np.diag(1 / np.diag(ds))
    idr = np.diag(1 / np.diag(dr))

    cgpts, residuals, relative = [], [], []
    for data in msr_list:
        cgpt = 4 * ids @ cs.T @ data @ cr @ idr / (nb_src * nb_rcv)
        res = float(np.linalg.norm(data - cs @ ds @ cgpt @ (cr @ dr).T))
        cgpts.append(cgpt)
        residuals.append(res)
        relative.append(res / float(np.linalg.norm(data)))

    return ReconstructionResult(
        cgpt=cgpts,
        residual=residuals,
        relative_residual=relative,
        source_matrix=cs @ ds,
        receiver_matrix=cr @ dr,
    )

plot

plot(ax=None, **kwargs)

Plot the inclusions and the acquisition system.

Parameters:

Name	Type	Description	Default
ax	Axes	Axes to draw on. A new figure is created if omitted.	None
**kwargs		Forwarded to the inclusion plot calls.	{}

Returns:

Type	Description
Axes	The axes containing the plot.

Source code in src/sies/pde/conductivity.py

def plot(self, ax=None, **kwargs):
    """Plot the inclusions and the acquisition system.

    Parameters
    ----------
    ax : matplotlib.axes.Axes, optional
        Axes to draw on. A new figure is created if omitted.
    **kwargs
        Forwarded to the inclusion plot calls.

    Returns
    -------
    matplotlib.axes.Axes
        The axes containing the plot.
    """
    import matplotlib.pyplot as plt

    if ax is None:
        _, ax = plt.subplots()
    for incl in self.inclusions:
        incl.plot(ax=ax, **kwargs)
    self.cfg.plot(ax=ax)
    return ax