Asymptotics (CGPT)

sies.asymptotics.cgpt

Numerical computation of contracted generalized polarization tensors.

The CGPT matrix of order k of inclusions \(D_1, \dots, D_L\) with contrasts \(\lambda_l\) is built from the boundary densities

$$ \phi_m = (\lambda I - K_D^)^{-1} \left[\frac{\partial \mathrm{Re}/\mathrm{Im}(z^m)}{\partial\nu}\right], $$ see Ammari & Kang, Polarization and Moment Tensors* (2007), chapter 4.

contrast

contrast(cnd, pmtt=None, freq=0.0)

Compute the contrast constants \(\lambda\) of inclusions.

For conductivity k, permittivity eps and frequency f, the contrast is

\[\lambda = \frac{k + 2\pi i f \epsilon + 1} {2(k + 2\pi i f \epsilon - 1)}.\]

Parameters:

Name	Type	Description	Default
cnd	array_like	Conductivity of each inclusion; must be positive and different from one (the background conductivity).	required
pmtt	array_like	Permittivity of each inclusion. Defaults to zero.	None
freq	float	Working frequency; must be nonnegative.	0.0

Returns:

Type	Description
ndarray	Complex contrast constant of each inclusion.

Notes

The \(2\pi\) factor comes from the Fourier transform convention \(\hat{f}(\omega) = \int f(x) e^{-2\pi i x \omega} dx\) used throughout the library (compatible with the FFT); published papers use the convention without the factor.

Source code in src/sies/asymptotics/cgpt.py

def contrast(
    cnd: NDArray,
    pmtt: NDArray | None = None,
    freq: float = 0.0,
) -> NDArray:
    r"""Compute the contrast constants $\lambda$ of inclusions.

    For conductivity ``k``, permittivity ``eps`` and frequency ``f``, the
    contrast is

    $$\lambda = \frac{k + 2\pi i f \epsilon + 1} {2(k + 2\pi i f \epsilon - 1)}.$$

    Parameters
    ----------
    cnd : array_like
        Conductivity of each inclusion; must be positive and different
        from one (the background conductivity).
    pmtt : array_like, optional
        Permittivity of each inclusion. Defaults to zero.
    freq : float, default 0.0
        Working frequency; must be nonnegative.

    Returns
    -------
    ndarray
        Complex contrast constant of each inclusion.

    Notes
    -----
    The $2\pi$ factor comes from the Fourier transform convention
    $\hat{f}(\omega) = \int f(x) e^{-2\pi i x \omega} dx$ used
    throughout the library (compatible with the FFT); published papers
    use the convention without the factor.
    """
    cnd = np.atleast_1d(np.asarray(cnd, dtype=float))
    if pmtt is None:
        pmtt = np.zeros_like(cnd)
    pmtt = np.atleast_1d(np.asarray(pmtt, dtype=float))

    if freq < 0:
        raise ValueError("Frequency must be a nonnegative scalar.")
    if np.any(cnd == 1) or np.any(cnd < 0):
        raise ValueError("Conductivity must be positive and different from 1.")

    kappa = cnd + 2j * np.pi * pmtt * freq
    lambdas = (kappa + 1) / (2 * (kappa - 1))
    return np.real_if_close(lambdas)

make_block_matrix

make_block_matrix(inclusions)

Assemble the frequency-independent blocks of the system matrix.

The diagonal blocks are \(-K_{D_n}^*\) and the off-diagonal (m, n) blocks are \(-\partial_{\nu_m} S_{D_n}\). The contrast-dependent term \(\lambda_n I\) is added later by make_system_matrix, which makes multi-frequency computations cheap.

Parameters:

Name	Type	Description	Default
inclusions	list of C2Boundary	Mutually disjoint inclusions, all discretized with the same number of boundary points.	required

Returns:

Type	Description
(ndarray, shape(L * n, L * n))	The block matrix, where L is the number of inclusions and n the number of boundary points.

Source code in src/sies/asymptotics/cgpt.py

def make_block_matrix(inclusions: list[C2Boundary]) -> NDArray:
    r"""Assemble the frequency-independent blocks of the system matrix.

    The diagonal blocks are $-K_{D_n}^*$ and the off-diagonal
    ``(m, n)`` blocks are
    $-\partial_{\nu_m} S_{D_n}$. The contrast-dependent term
    $\lambda_n I$ is added later by `make_system_matrix`,
    which makes multi-frequency computations cheap.

    Parameters
    ----------
    inclusions : list of C2Boundary
        Mutually disjoint inclusions, all discretized with the same
        number of boundary points.

    Returns
    -------
    ndarray, shape (L * n, L * n)
        The block matrix, where ``L`` is the number of inclusions and
        ``n`` the number of boundary points.
    """
    _check_disjoint(inclusions)
    nb_incl = len(inclusions)
    nb_points = inclusions[0].nb_points
    blocks = np.empty((nb_incl * nb_points, nb_incl * nb_points))

    for m, image in enumerate(inclusions):
        rows = slice(m * nb_points, (m + 1) * nb_points)
        for n, domain in enumerate(inclusions):
            cols = slice(n * nb_points, (n + 1) * nb_points)
            if m == n:
                blocks[rows, cols] = -KStar(domain).matrix
            else:
                blocks[rows, cols] = -SingleLayerNormalDerivative.kernel_matrix(
                    domain.points, domain.sigma, image.points, image.normal
                )
    return blocks

make_system_matrix

make_system_matrix(block_matrix, lambdas)

Assemble the full system matrix \(\lambda I - K_D^*\).

Parameters:

Name	Type	Description	Default
block_matrix	ndarray	Frequency-independent blocks from make_block_matrix.	required
lambdas	array_like	Contrast constant of each inclusion.	required

Returns:

Type	Description
ndarray	The system matrix A of the linear system A phi = b.

Source code in src/sies/asymptotics/cgpt.py

def make_system_matrix(block_matrix: NDArray, lambdas: NDArray) -> NDArray:
    r"""Assemble the full system matrix $\lambda I - K_D^*$.

    Parameters
    ----------
    block_matrix : ndarray
        Frequency-independent blocks from `make_block_matrix`.
    lambdas : array_like
        Contrast constant of each inclusion.

    Returns
    -------
    ndarray
        The system matrix ``A`` of the linear system ``A phi = b``.
    """
    lambdas = np.atleast_1d(np.asarray(lambdas))
    nb_incl = len(lambdas)
    nb_points = block_matrix.shape[0] // nb_incl

    system = block_matrix.astype(np.result_type(block_matrix, lambdas), copy=True)
    for n, lam in enumerate(lambdas):
        idx = np.arange(n * nb_points, (n + 1) * nb_points)
        system[idx, idx] += lam
    return system

theoretical_cgpt

theoretical_cgpt(inclusions, lambdas, order, block_matrix=None)

Compute the CGPT matrix of inclusions up to a given order.

Entries follow the convention

$$ M^{cc}{mn} = \int{\partial D} \mathrm{Re}(z^n) \, (\lambda I - K_D^*)^{-1} \left[\partial_\nu \mathrm{Re}(z^m)\right] ds, $$ and similarly for the cs, sc and ss blocks, interleaved in a (2 * order, 2 * order) matrix.

Parameters:

Name	Type	Description	Default
inclusions	C2Boundary or list of C2Boundary	The inclusion(s).	required
lambdas	array_like	Contrast constant of each inclusion.	required
order	int	Maximum CGPT order.	required
block_matrix	ndarray	Precomputed blocks from make_block_matrix, to be reused across frequencies.	None

Returns:

Type	Description
(ndarray, shape(2 * order, 2 * order))	The CGPT matrix (complex if any contrast is complex).

Source code in src/sies/asymptotics/cgpt.py

def theoretical_cgpt(
    inclusions: list[C2Boundary] | C2Boundary,
    lambdas: NDArray,
    order: int,
    block_matrix: NDArray | None = None,
) -> NDArray:
    r"""Compute the CGPT matrix of inclusions up to a given order.

    Entries follow the convention

    $$
    M^{cc}_{mn} = \int_{\partial D} \mathrm{Re}(z^n) \,
    (\lambda I - K_D^*)^{-1} \left[\partial_\nu \mathrm{Re}(z^m)\right] ds,
    $$
    and similarly for the ``cs``, ``sc`` and ``ss`` blocks, interleaved
    in a ``(2 * order, 2 * order)`` matrix.

    Parameters
    ----------
    inclusions : C2Boundary or list of C2Boundary
        The inclusion(s).
    lambdas : array_like
        Contrast constant of each inclusion.
    order : int
        Maximum CGPT order.
    block_matrix : ndarray, optional
        Precomputed blocks from `make_block_matrix`, to be reused
        across frequencies.

    Returns
    -------
    ndarray, shape (2 * order, 2 * order)
        The CGPT matrix (complex if any contrast is complex).
    """
    if isinstance(inclusions, C2Boundary):
        inclusions = [inclusions]
    lambdas = np.atleast_1d(np.asarray(lambdas))
    if len(lambdas) < len(inclusions):
        raise ValueError("A contrast value must be specified for each inclusion.")

    if block_matrix is None:
        block_matrix = make_block_matrix(inclusions)

    nb_incl = len(inclusions)
    nb_points = inclusions[0].nb_points
    system = make_system_matrix(block_matrix, lambdas)

    # When lambda is close to 1/2 (infinite contrast), enforce that the
    # densities have zero mean by augmenting the system.
    augmented = bool(np.min(np.abs(lambdas - 0.5)) < _HALF_CONTRAST_TOL)
    if augmented:
        constraints = np.kron(np.eye(nb_incl), np.ones((1, nb_points)))
        system = np.vstack([system, constraints])

    # Right-hand sides nu . grad(z^m) for all orders m, stacked by inclusion.
    rhs = np.empty((nb_incl * nb_points, order), dtype=complex)
    for i, incl in enumerate(inclusions):
        rows = slice(i * nb_points, (i + 1) * nb_points)
        for m in range(1, order + 1):
            grad = m * incl.cpoints ** (m - 1)
            rhs[rows, m - 1] = (incl.normal[0] + 1j * incl.normal[1]) * grad
    if augmented:
        rhs = np.vstack([rhs, np.zeros((nb_incl, order))])

    phi_re = _solve_densities(system, rhs.real, augmented)
    phi_im = _solve_densities(system, rhs.imag, augmented)

    is_complex = np.iscomplexobj(phi_re) or np.iscomplexobj(phi_im)
    dtype = complex if is_complex else float
    cc = np.zeros((order, order), dtype=dtype)
    cs = np.zeros_like(cc)
    sc = np.zeros_like(cc)
    ss = np.zeros_like(cc)

    for i, incl in enumerate(inclusions):
        rows = slice(i * nb_points, (i + 1) * nb_points)
        # moments[n - 1] = z^n * sigma for n = 1..order, shape (order, nb_points)
        powers = incl.cpoints[np.newaxis, :] ** np.arange(1, order + 1)[:, np.newaxis]
        moments = powers * incl.sigma
        cc += moments.real @ phi_re[rows]
        cs += moments.imag @ phi_re[rows]
        sc += moments.real @ phi_im[rows]
        ss += moments.imag @ phi_im[rows]

    # The (m, n) entry pairs the m-th density with the n-th moment.
    return join_cgpt(cc.T, cs.T, sc.T, ss.T)

sies.asymptotics.exact

Closed-form CGPT matrices of disks and ellipses.

These exact formulas (M. Lim's PhD thesis, also reported in Ammari & Kang 2007) serve as reference values to validate the boundary integral computation of theoretical_cgpt.

disk_cgpt

disk_cgpt(order, radius, cnd)

Exact CGPT matrix of a disk centered at the origin.

Parameters:

Name	Type	Description	Default
order	int	Maximum CGPT order.	required
radius	float	Radius of the disk.	required
cnd	float	Conductivity of the disk (background is one).	required

Returns:

Type	Description
(ndarray, shape(2 * order, 2 * order))	The diagonal CGPT matrix of the disk.

Source code in src/sies/asymptotics/exact.py

def disk_cgpt(order: int, radius: float, cnd: float) -> NDArray:
    """Exact CGPT matrix of a disk centered at the origin.

    Parameters
    ----------
    order : int
        Maximum CGPT order.
    radius : float
        Radius of the disk.
    cnd : float
        Conductivity of the disk (background is one).

    Returns
    -------
    ndarray, shape (2 * order, 2 * order)
        The diagonal CGPT matrix of the disk.
    """
    n = np.arange(1, order + 1)
    diag = np.pi * n * radius ** (2 * n) * 2 * (cnd - 1) / (cnd + 1)
    return np.diag(np.repeat(diag, 2))

ellipse_cgpt

ellipse_cgpt(order, axis_a, axis_b, cnd)

Exact CGPT matrix of an ellipse centered at the origin.

Parameters:

Name	Type	Description	Default
order	int	Maximum CGPT order.	required
axis_a	float	Semi-major axis length.	required
axis_b	float	Semi-minor axis length.	required
cnd	float	Conductivity of the ellipse (background is one). Complex values are not supported by the closed form.	required

Returns:

Type	Description
(ndarray, shape(2 * order, 2 * order))	The CGPT matrix of the ellipse.

Source code in src/sies/asymptotics/exact.py

def ellipse_cgpt(order: int, axis_a: float, axis_b: float, cnd: float) -> NDArray:
    """Exact CGPT matrix of an ellipse centered at the origin.

    Parameters
    ----------
    order : int
        Maximum CGPT order.
    axis_a : float
        Semi-major axis length.
    axis_b : float
        Semi-minor axis length.
    cnd : float
        Conductivity of the ellipse (background is one). Complex values
        are not supported by the closed form.

    Returns
    -------
    ndarray, shape (2 * order, 2 * order)
        The CGPT matrix of the ellipse.
    """
    exact = np.zeros((2 * order, 2 * order))
    for m in range(1, order + 1):
        for n in range(m, order + 1):
            block = _ellipse_block(m, n, axis_a, axis_b, cnd)
            exact[2 * m - 2 : 2 * m, 2 * n - 2 : 2 * n] = block.real
    # Fill the lower triangle by symmetry of the CGPT.
    upper = np.triu(np.ones_like(exact), 1)
    exact = exact + (exact * upper).T
    return exact * (cnd - 1)

sies.asymptotics.transforms

Transformation rules of CGPT matrices under rigid motions and scaling.

A CGPT matrix is stored as a (2k, 2k) array interleaving the four cc, cs, sc and ss blocks. The complex CGPT pair (N1, N2) diagonalizes the action of translations, rotations and dilations, which makes these transformations explicit; see Ammari et al., Target identification using dictionary matching of generalized polarization tensors (2014).

split_cgpt

split_cgpt(cgpt)

Split an interleaved CGPT matrix into its four blocks.

Parameters:

Name	Type	Description	Default
cgpt	ndarray, shape (2k, 2k)	Interleaved CGPT matrix.	required

Returns:

Type	Description
cc, cs, sc, ss : ndarray, shape (k, k)	The four blocks of the CGPT.

Source code in src/sies/asymptotics/transforms.py

def split_cgpt(cgpt: NDArray) -> tuple[NDArray, NDArray, NDArray, NDArray]:
    """Split an interleaved CGPT matrix into its four blocks.

    Parameters
    ----------
    cgpt : ndarray, shape (2k, 2k)
        Interleaved CGPT matrix.

    Returns
    -------
    cc, cs, sc, ss : ndarray, shape (k, k)
        The four blocks of the CGPT.
    """
    return cgpt[0::2, 0::2], cgpt[0::2, 1::2], cgpt[1::2, 0::2], cgpt[1::2, 1::2]

join_cgpt

join_cgpt(cc, cs, sc, ss)

Interleave the four CGPT blocks into a single matrix.

Parameters:

Name	Type	Description	Default
cc	(ndarray, shape(k, k))	The four blocks of the CGPT.	required
cs	(ndarray, shape(k, k))	The four blocks of the CGPT.	required
sc	(ndarray, shape(k, k))	The four blocks of the CGPT.	required
ss	(ndarray, shape(k, k))	The four blocks of the CGPT.	required

Returns:

Type	Description
ndarray, shape (2k, 2k)	The interleaved CGPT matrix.

Source code in src/sies/asymptotics/transforms.py

def join_cgpt(cc: NDArray, cs: NDArray, sc: NDArray, ss: NDArray) -> NDArray:
    """Interleave the four CGPT blocks into a single matrix.

    Parameters
    ----------
    cc, cs, sc, ss : ndarray, shape (k, k)
        The four blocks of the CGPT.

    Returns
    -------
    ndarray, shape (2k, 2k)
        The interleaved CGPT matrix.
    """
    order = cc.shape[0]
    cgpt = np.zeros((2 * order, 2 * order), dtype=np.result_type(cc, cs, sc, ss))
    cgpt[0::2, 0::2] = cc
    cgpt[0::2, 1::2] = cs
    cgpt[1::2, 0::2] = sc
    cgpt[1::2, 1::2] = ss
    return cgpt

cgpt_to_complex

cgpt_to_complex(cgpt)

Build the complex CGPT pair from an interleaved CGPT matrix.

Parameters:

Name	Type	Description	Default
cgpt	ndarray, shape (2k, 2k)	Interleaved CGPT matrix.	required

Returns:

Name	Type	Description
n1	(ndarray, shape(k, k))	First complex CGPT; symmetric.
n2	(ndarray, shape(k, k))	Second complex CGPT; hermitian.

Source code in src/sies/asymptotics/transforms.py

def cgpt_to_complex(cgpt: NDArray) -> tuple[NDArray, NDArray]:
    """Build the complex CGPT pair from an interleaved CGPT matrix.

    Parameters
    ----------
    cgpt : ndarray, shape (2k, 2k)
        Interleaved CGPT matrix.

    Returns
    -------
    n1 : ndarray, shape (k, k)
        First complex CGPT; symmetric.
    n2 : ndarray, shape (k, k)
        Second complex CGPT; hermitian.
    """
    cc, cs, sc, ss = split_cgpt(cgpt)
    n1 = cc - ss + 1j * (cs + sc)
    n2 = cc + ss + 1j * (cs - sc)
    return n1, n2

complex_to_cgpt

complex_to_cgpt(n1, n2)

Convert the complex CGPT pair back to an interleaved CGPT matrix.

Only valid for real contrasts (real-valued CGPT).

Parameters:

Name	Type	Description	Default
n1	(ndarray, shape(k, k))	First complex CGPT.	required
n2	(ndarray, shape(k, k))	Second complex CGPT.	required

Returns:

Type	Description
ndarray, shape (2k, 2k)	The interleaved CGPT matrix.

Source code in src/sies/asymptotics/transforms.py

def complex_to_cgpt(n1: NDArray, n2: NDArray) -> NDArray:
    """Convert the complex CGPT pair back to an interleaved CGPT matrix.

    Only valid for real contrasts (real-valued CGPT).

    Parameters
    ----------
    n1 : ndarray, shape (k, k)
        First complex CGPT.
    n2 : ndarray, shape (k, k)
        Second complex CGPT.

    Returns
    -------
    ndarray, shape (2k, 2k)
        The interleaved CGPT matrix.
    """
    cc = (n1 + n2).real / 2
    cs = (n1 + n2).imag / 2
    sc = (n1 - n2).imag / 2
    ss = (n2 - n1).real / 2
    return join_cgpt(cc, cs, sc, ss)

transform_ccgpt

transform_ccgpt(n1, n2, translation, scaling=1.0, rotation=0.0)

Apply rotation, then scaling, then translation to a complex CGPT pair.

Parameters:

Name	Type	Description	Default
n1	(ndarray, shape(k, k))	Complex CGPT pair of the original shape.	required
n2	(ndarray, shape(k, k))	Complex CGPT pair of the original shape.	required
translation	complex or array_like	Translation, as a complex number tx + i ty or a 2-vector.	required
scaling	float	Scaling factor.	1.0
rotation	float	Rotation angle in radians.	0.0

Returns:

Type	Description
z1, z2 : ndarray, shape (k, k)	Complex CGPT pair of the transformed shape.

Source code in src/sies/asymptotics/transforms.py

def transform_ccgpt(
    n1: NDArray,
    n2: NDArray,
    translation: complex,
    scaling: float = 1.0,
    rotation: float = 0.0,
) -> tuple[NDArray, NDArray]:
    """Apply rotation, then scaling, then translation to a complex CGPT pair.

    Parameters
    ----------
    n1, n2 : ndarray, shape (k, k)
        Complex CGPT pair of the original shape.
    translation : complex or array_like
        Translation, as a complex number ``tx + i ty`` or a 2-vector.
    scaling : float, default 1.0
        Scaling factor.
    rotation : float, default 0.0
        Rotation angle in radians.

    Returns
    -------
    z1, z2 : ndarray, shape (k, k)
        Complex CGPT pair of the transformed shape.
    """
    order = n1.shape[0]
    t0 = _as_complex(translation)
    ct = _translation_matrix(t0, order, inverse=False)
    gy = np.diag((scaling * np.exp(1j * rotation)) ** np.arange(1, order + 1))

    z1 = ct @ gy @ n1 @ gy @ ct.T
    z2 = ct.conj() @ gy.conj() @ n2 @ gy @ ct.T
    return z1, z2

transform_ccgpt_inverse

transform_ccgpt_inverse(n1, n2, translation, scaling=1.0, rotation=0.0)

Apply the inverse transformation to a complex CGPT pair.

Undo a translation, then a scaling, then a rotation (the inverse of transform_ccgpt with the same arguments).

Parameters:

Name	Type	Description	Default
n1	(ndarray, shape(k, k))	Complex CGPT pair of the transformed shape.	required
n2	(ndarray, shape(k, k))	Complex CGPT pair of the transformed shape.	required
translation	complex or array_like	Translation to undo.	required
scaling	float	Scaling factor to undo.	1.0
rotation	float	Rotation angle to undo, in radians.	0.0

Returns:

Type	Description
z1, z2 : ndarray, shape (k, k)	Complex CGPT pair with the transformation removed.

Source code in src/sies/asymptotics/transforms.py

def transform_ccgpt_inverse(
    n1: NDArray,
    n2: NDArray,
    translation: complex,
    scaling: float = 1.0,
    rotation: float = 0.0,
) -> tuple[NDArray, NDArray]:
    """Apply the inverse transformation to a complex CGPT pair.

    Undo a translation, then a scaling, then a rotation (the inverse of
    `transform_ccgpt` with the same arguments).

    Parameters
    ----------
    n1, n2 : ndarray, shape (k, k)
        Complex CGPT pair of the transformed shape.
    translation : complex or array_like
        Translation to undo.
    scaling : float, default 1.0
        Scaling factor to undo.
    rotation : float, default 0.0
        Rotation angle to undo, in radians.

    Returns
    -------
    z1, z2 : ndarray, shape (k, k)
        Complex CGPT pair with the transformation removed.
    """
    order = n1.shape[0]
    t0 = _as_complex(translation)
    ct = _translation_matrix(t0, order, inverse=True)
    gy = np.diag((np.exp(-1j * rotation) / scaling) ** np.arange(1, order + 1))

    z1 = gy @ ct @ n1 @ ct.T @ gy
    z2 = gy.conj() @ ct.conj() @ n2 @ ct.T @ gy
    return z1, z2

transform_cgpt

transform_cgpt(cgpt, translation, scaling=1.0, rotation=0.0)

Apply rotation, then scaling, then translation to a CGPT matrix.

Parameters:

Name	Type	Description	Default
cgpt	ndarray, shape (2k, 2k)	Interleaved CGPT matrix of the original shape.	required
translation	complex or array_like	Translation, as a complex number tx + i ty or a 2-vector.	required
scaling	float	Scaling factor.	1.0
rotation	float	Rotation angle in radians.	0.0

Returns:

Type	Description
ndarray, shape (2k, 2k)	CGPT matrix of the transformed shape.

Source code in src/sies/asymptotics/transforms.py

def transform_cgpt(
    cgpt: NDArray,
    translation: complex,
    scaling: float = 1.0,
    rotation: float = 0.0,
) -> NDArray:
    """Apply rotation, then scaling, then translation to a CGPT matrix.

    Parameters
    ----------
    cgpt : ndarray, shape (2k, 2k)
        Interleaved CGPT matrix of the original shape.
    translation : complex or array_like
        Translation, as a complex number ``tx + i ty`` or a 2-vector.
    scaling : float, default 1.0
        Scaling factor.
    rotation : float, default 0.0
        Rotation angle in radians.

    Returns
    -------
    ndarray, shape (2k, 2k)
        CGPT matrix of the transformed shape.
    """
    n1, n2 = cgpt_to_complex(cgpt)
    z1, z2 = transform_ccgpt(n1, n2, translation, scaling, rotation)
    return complex_to_cgpt(z1, z2)