__all__ = [
"zoeppritz_scattering",
"zoeppritz_element",
"zoeppritz_pp",
"approx_zoeppritz_pp",
"akirichards",
"fatti",
"ps",
"AVOLinearModelling",
]
import logging
from typing import List, Optional, Tuple, Union
import numpy as np
import numpy.typing as npt
from numpy import cos, sin, tan
from pylops import LinearOperator
from pylops.utils._internal import _value_or_sized_to_tuple
from pylops.utils.backend import get_array_module
from pylops.utils.decorators import reshaped
from pylops.utils.typing import DTypeLike, NDArray
logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.WARNING)
[docs]def zoeppritz_scattering(
vp1: float,
vs1: float,
rho1: float,
vp0: float,
vs0: float,
rho0: float,
theta1: Union[float, npt.ArrayLike],
) -> NDArray:
r"""Zoeppritz solution.
Calculates the angle dependent p-wave reflectivity of an
interface between two media for a set of incident angles.
Parameters
----------
vp1 : :obj:`float`
P-wave velocity of the upper medium
vs1 : :obj:`float`
S-wave velocity of the upper medium
rho1 : :obj:`float`
Density of the upper medium
vp0 : :obj:`float`
P-wave velocity of the lower medium
vs0 : :obj:`float`
S-wave velocity of the lower medium
rho0 : :obj:`float`
Density of the lower medium
theta1 : :obj:`np.ndarray` or :obj:`float`
Incident angles in degrees
Returns
-------
zoep : :obj:`np.ndarray`
:math:`4 \times 4` matrix representing the scattering matrix for the
incident angle ``theta1``
See also
--------
zoeppritz_element : Single reflectivity element of Zoeppritz solution
zoeppritz_pp : PP reflectivity element of Zoeppritz solution
"""
ncp = get_array_module(theta1)
# Create theta1 array of angles in radiants
if isinstance(theta1, (int, float)):
theta1 = ncp.array(
[
float(theta1),
]
)
elif isinstance(theta1, (list, tuple)):
theta1 = ncp.array(theta1)
theta1 = ncp.radians(theta1)
# Set the ray parameter p
p = sin(theta1) / vp1
# Calculate reflection & transmission angles for Zoeppritz
theta2 = ncp.arcsin(p * vp0) # Trans. angle of P-wave
phi1 = ncp.arcsin(p * vs1) # Refl. angle of converted S-wave
phi2 = ncp.arcsin(p * vs0) # Trans. angle of converted S-wave
# Matrix form of Zoeppritz equation
M = ncp.array(
[
[-sin(theta1), -cos(phi1), sin(theta2), cos(phi2)],
[cos(theta1), -sin(phi1), cos(theta2), -sin(phi2)],
[
2 * rho1 * vs1 * sin(phi1) * cos(theta1),
rho1 * vs1 * (1 - 2 * sin(phi1) ** 2),
2 * rho0 * vs0 * sin(phi2) * cos(theta2),
rho0 * vs0 * (1 - 2 * sin(phi2) ** 2),
],
[
-rho1 * vp1 * (1 - 2 * sin(phi1) ** 2),
rho1 * vs1 * sin(2 * phi1),
rho0 * vp0 * (1 - 2 * sin(phi2) ** 2),
-rho0 * vs0 * sin(2 * phi2),
],
],
dtype="float",
)
N = ncp.array(
[
[sin(theta1), cos(phi1), -sin(theta2), -cos(phi2)],
[cos(theta1), -sin(phi1), cos(theta2), -sin(phi2)],
[
2 * rho1 * vs1 * sin(phi1) * cos(theta1),
rho1 * vs1 * (1 - 2 * sin(phi1) ** 2),
2 * rho0 * vs0 * sin(phi2) * cos(theta2),
rho0 * vs0 * (1 - 2 * sin(phi2) ** 2),
],
[
rho1 * vp1 * (1 - 2 * sin(phi1) ** 2),
-rho1 * vs1 * sin(2 * phi1),
-rho0 * vp0 * (1 - 2 * sin(phi2) ** 2),
rho0 * vs0 * sin(2 * phi2),
],
],
dtype="float",
)
# Create Zoeppritz coefficient for all angles
zoep = ncp.zeros((4, 4, M.shape[-1]))
for i in range(M.shape[-1]):
Mi = M[..., i]
Ni = N[..., i]
dt = ncp.dot(ncp.linalg.inv(Mi), Ni)
zoep[..., i] = dt
return zoep
[docs]def zoeppritz_element(
vp1: float,
vs1: float,
rho1: float,
vp0: float,
vs0: float,
rho0: float,
theta1: Union[float, NDArray],
element: str = "PdPu",
) -> NDArray:
"""Single element of Zoeppritz solution.
Simple wrapper to :py:class:`pylops.avo.avo.scattering_matrix`,
returning any mode reflection coefficient from the Zoeppritz
scattering matrix for specific combination of incident
and reflected wave and a set of incident angles
Parameters
----------
vp1 : :obj:`float`
P-wave velocity of the upper medium
vs1 : :obj:`float`
S-wave velocity of the upper medium
rho1 : :obj:`float`
Density of the upper medium
vp0 : :obj:`float`
P-wave velocity of the lower medium
vs0 : :obj:`float`
S-wave velocity of the lower medium
rho0 : :obj:`float`
Density of the lower medium
theta1 : :obj:`np.ndarray` or :obj:`float`
Incident angles in degrees
element : :obj:`str`, optional
Specific choice of incident and reflected wave combining
any two of the following strings: ``Pd`` P-wave downgoing,
``Sd`` S-wave downgoing, ``Pu`` P-wave upgoing,
``Su`` S-wave upgoing (e.g., ``PdPu``)
Returns
-------
refl : :obj:`np.ndarray`
reflectivity values for all input angles for specific combination
of incident and reflected wave.
See also
--------
zoeppritz_scattering : Zoeppritz solution
zoeppritz_pp : PP reflectivity element of Zoeppritz solution
"""
elements = np.array(
[
["PdPu", "SdPu", "PuPu", "SuPu"],
["PdSu", "SdSu", "PuSu", "SuSu"],
["PdPd", "SdPd", "PuPd", "SuPd"],
["PdSd", "SdSd", "PuSd", "SuSd"],
]
)
refl = zoeppritz_scattering(vp1, vs1, rho1, vp0, vs0, rho0, theta1)
element = np.where(elements == element)
return np.squeeze(refl[element])
[docs]def zoeppritz_pp(
vp1: float,
vs1: float,
rho1: float,
vp0: float,
vs0: float,
rho0: float,
theta1: Union[float, NDArray],
) -> NDArray:
"""PP reflection coefficient from the Zoeppritz scattering matrix.
Simple wrapper to :py:class:`pylops.avo.avo.scattering_matrix`,
returning the PP reflection coefficient from the Zoeppritz
scattering matrix for a set of incident angles
Parameters
----------
vp1 : :obj:`float`
P-wave velocity of the upper medium
vs1 : :obj:`float`
S-wave velocity of the upper medium
rho1 : :obj:`float`
Density of the upper medium
vp0 : :obj:`float`
P-wave velocity of the lower medium
vs0 : :obj:`float`
S-wave velocity of the lower medium
rho0 : :obj:`float`
Density of the lower medium
theta1 : :obj:`np.ndarray` or :obj:`float`
Incident angles in degrees
Returns
-------
PPrefl : :obj:`np.ndarray`
PP reflectivity values for all input angles.
See also
--------
zoeppritz_scattering : Zoeppritz solution
zoeppritz_element : Single reflectivity element of Zoeppritz solution
"""
PPrefl = zoeppritz_element(vp1, vs1, rho1, vp0, vs0, rho0, theta1, "PdPu")
return PPrefl
[docs]def approx_zoeppritz_pp(
vp1: Union[List, Tuple, npt.ArrayLike],
vs1: Union[List, Tuple, npt.ArrayLike],
rho1: Union[List, Tuple, npt.ArrayLike],
vp0: Union[List, Tuple, npt.ArrayLike],
vs0: Union[List, Tuple, npt.ArrayLike],
rho0: Union[List, Tuple, npt.ArrayLike],
theta1: Union[float, NDArray],
) -> NDArray:
"""PP reflection coefficient from the approximate Zoeppritz equation.
Approximate calculation of PP reflection from the Zoeppritz
scattering matrix for a set of incident angles [1]_.
.. [1] Dvorkin et al. Seismic Reflections of Rock Properties.
Cambridge. 2014.
Parameters
----------
vp1 : :obj:`np.ndarray` or :obj:`list` or :obj:`tuple`
P-wave velocity of the upper medium
vs1 : :obj:`np.ndarray` or :obj:`list` or :obj:`tuple`
S-wave velocity of the upper medium
rho1 : :obj:`np.ndarray` or :obj:`list` or :obj:`tuple`
Density of the upper medium
vp0 : :obj:`np.ndarray` or :obj:`list` or :obj:`tuple`
P-wave velocity of the lower medium
vs0 : :obj:`np.ndarray` or :obj:`list` or :obj:`tuple`
S-wave velocity of the lower medium
rho0 : :obj:`np.ndarray` or :obj:`list` or :obj:`tuple`
Density of the lower medium
theta1 : :obj:`np.ndarray` or :obj:`float`
Incident angles in degrees
Returns
-------
PPrefl : :obj:`np.ndarray`
PP reflectivity values for all input angles.
See also
--------
zoeppritz_scattering : Zoeppritz solution
zoeppritz_element : Single reflectivity element of Zoeppritz solution
zoeppritz_pp : PP reflectivity element of Zoeppritz solution
"""
ncp = get_array_module(theta1)
vp1, vs1, rho1 = ncp.array(vp1), ncp.array(vs1), ncp.array(rho1)
vp0, vs0, rho0 = ncp.array(vp0), ncp.array(vs0), ncp.array(rho0)
# Incident P
theta1 = theta1[:, np.newaxis] if vp1.size > 1 else theta1
theta1 = ncp.deg2rad(theta1)
# Ray parameter and reflected P
p = ncp.sin(theta1) / vp1
theta0 = ncp.arcsin(p * vp0)
# Reflected S
phi1 = ncp.arcsin(p * vs1)
# Transmitted S
phi0 = ncp.arcsin(p * vs0)
# Coefficients
a = rho0 * (1 - 2 * np.sin(phi0) ** 2.0) - rho1 * (1 - 2 * np.sin(phi1) ** 2.0)
b = rho0 * (1 - 2 * np.sin(phi0) ** 2.0) + 2 * rho1 * np.sin(phi1) ** 2.0
c = rho1 * (1 - 2 * np.sin(phi1) ** 2.0) + 2 * rho0 * np.sin(phi0) ** 2.0
d = 2 * (rho0 * vs0**2 - rho1 * vs1**2)
E = (b * np.cos(theta1) / vp1) + (c * np.cos(theta0) / vp0)
F = (b * np.cos(phi1) / vs1) + (c * np.cos(phi0) / vs0)
G = a - d * np.cos(theta1) / vp1 * np.cos(phi0) / vs0
H = a - d * np.cos(theta0) / vp0 * np.cos(phi1) / vs1
D = E * F + G * H * p**2
rpp = (1 / D) * (
F * (b * (ncp.cos(theta1) / vp1) - c * (ncp.cos(theta0) / vp0))
- H * p**2 * (a + d * (ncp.cos(theta1) / vp1) * (ncp.cos(phi0) / vs0))
)
return rpp
[docs]def akirichards(
theta: npt.ArrayLike,
vsvp: Union[float, NDArray],
n: int = 1,
) -> Tuple[NDArray, NDArray, NDArray]:
r"""Three terms Aki-Richards approximation.
Computes the coefficients of the of three terms Aki-Richards approximation
for a set of angles and a constant or variable VS/VP ratio.
Parameters
----------
theta : :obj:`np.ndarray`
Incident angles in degrees
vsvp : :obj:`np.ndarray` or :obj:`float`
:math:`V_S/V_P` ratio
n : :obj:`int`, optional
Number of samples (if ``vsvp`` is a scalar)
Returns
-------
G1 : :obj:`np.ndarray`
First coefficient of three terms Aki-Richards approximation
:math:`[n_\theta \times n_\text{vsvp}]`
G2 : :obj:`np.ndarray`
Second coefficient of three terms Aki-Richards approximation
:math:`[n_\theta \times n_\text{vsvp}]`
G3 : :obj:`np.ndarray`
Third coefficient of three terms Aki-Richards approximation
:math:`[n_\theta \times n_\text{vsvp}]`
Notes
-----
The three terms Aki-Richards approximation [1]_, [2]_, is used to compute the
reflection coefficient as linear combination of contrasts in
:math:`V_P`, :math:`V_S`, and :math:`\rho.` More specifically:
.. math::
R(\theta) = G_1(\theta) \frac{\Delta V_P}{\overline{V}_P} + G_2(\theta)
\frac{\Delta V_S}{\overline{V}_S} + G_3(\theta)
\frac{\Delta \rho}{\overline{\rho}}
where
.. math::
\begin{align}
G_1(\theta) &= \frac{1}{2 \cos^2 \theta},\\
G_2(\theta) &= -4 (V_S/V_P)^2 \sin^2 \theta,\\
G_3(\theta) &= 0.5 - 2 (V_S/V_P)^2 \sin^2 \theta,\\
\frac{\Delta V_P}{\overline{V}_P} &= 2 \frac{V_{P,2}-V_{P,1}}{V_{P,2}+V_{P,1}},\\
\frac{\Delta V_S}{\overline{V}_S} &= 2 \frac{V_{S,2}-V_{S,1}}{V_{S,2}+V_{S,1}}, \\
\frac{\Delta \rho}{\overline{\rho}} &= 2 \frac{\rho_2-\rho_1}{\rho_2+\rho_1}.
\end{align}
.. [1] https://wiki.seg.org/wiki/AVO_equations
.. [2] Aki, K., and Richards, P. G. (2002). Quantitative Seismology (2nd ed.). University Science Books.
"""
ncp = get_array_module(theta)
theta = ncp.deg2rad(theta)
vsvp = vsvp * ncp.ones(n) if not isinstance(vsvp, ncp.ndarray) else vsvp
theta = theta[:, np.newaxis] if vsvp.size > 1 else theta
vsvp = vsvp[:, np.newaxis].T if vsvp.size > 1 else vsvp
G1 = 1.0 / (2.0 * cos(theta) ** 2) + 0 * vsvp
G2 = -4.0 * vsvp**2 * np.sin(theta) ** 2
G3 = 0.5 - 2.0 * vsvp**2 * sin(theta) ** 2
return G1, G2, G3
[docs]def fatti(
theta: npt.ArrayLike,
vsvp: Union[float, NDArray],
n: int = 1,
) -> Tuple[NDArray, NDArray, NDArray]:
r"""Three terms Fatti approximation.
Computes the coefficients of the three terms Fatti approximation
for a set of angles and a constant or variable VS/VP ratio.
Parameters
----------
theta : :obj:`np.ndarray`
Incident angles in degrees
vsvp : :obj:`np.ndarray` or :obj:`float`
:math:`V_S/V_P` ratio
n : :obj:`int`, optional
Number of samples (if ``vsvp`` is a scalar)
Returns
-------
G1 : :obj:`np.ndarray`
First coefficient of three terms Smith-Gidlow approximation
:math:`[n_{\theta} \times n_\text{vsvp}]`
G2 : :obj:`np.ndarray`
Second coefficient of three terms Smith-Gidlow approximation
:math:`[n_{\theta} \times n_\text{vsvp}]`
G3 : :obj:`np.ndarray`
Third coefficient of three terms Smith-Gidlow approximation
:math:`[n_{\theta} \times n_\text{vsvp}]`
Notes
-----
The three terms Fatti approximation [1]_, [2]_, is used to compute the reflection
coefficient as linear combination of contrasts in :math:`\text{AI},`
:math:`\text{SI}`, and :math:`\rho.` More specifically:
.. math::
R(\theta) = G_1(\theta) \frac{\Delta \text{AI}}{\bar{\text{AI}}} + G_2(\theta)
\frac{\Delta \text{SI}}{\overline{\text{SI}}} +
G_3(\theta) \frac{\Delta \rho}{\overline{\rho}}
where
.. math::
\begin{align}
G_1(\theta) &= 0.5 (1 + \tan^2 \theta),\\
G_2(\theta) &= -4 (V_S/V_P)^2 \sin^2 \theta,\\
G_3(\theta) &= 0.5 \left(4 (V_S/V_P)^2 \sin^2 \theta - \tan^2 \theta\right),\\
\frac{\Delta \text{AI}}{\overline{\text{AI}}} &= 2 \frac{\text{AI}_2-\text{AI}_1}{\text{AI}_2+\text{AI}_1},\\
\frac{\Delta \text{SI}}{\overline{\text{SI}}} &= 2 \frac{\text{SI}_2-\text{SI}_1}{\text{SI}_2+\text{SI}_1},\\
\frac{\Delta \rho}{\overline{\rho}} &= 2 \frac{\rho_2-\rho_1}{\rho_2+\rho_1}.
\end{align}
.. [1] https://www.subsurfwiki.org/wiki/Fatti_equation
.. [2] Jan L. Fatti, George C. Smith, Peter J. Vail, Peter J. Strauss, and Philip R. Levitt, (1994), "Detection of gas in sandstone reservoirs using AVO analysis: A 3-D seismic case history using the Geostack technique," Geophysics 59: 1362-1376.
"""
ncp = get_array_module(theta)
theta = ncp.deg2rad(theta)
vsvp = vsvp * ncp.ones(n) if not isinstance(vsvp, ncp.ndarray) else vsvp
theta = theta[:, np.newaxis] if vsvp.size > 1 else theta
vsvp = vsvp[:, np.newaxis].T if vsvp.size > 1 else vsvp
G1 = 0.5 * (1 + np.tan(theta) ** 2) + 0 * vsvp
G2 = -4 * vsvp**2 * np.sin(theta) ** 2
G3 = 0.5 * (4 * vsvp**2 * np.sin(theta) ** 2 - tan(theta) ** 2)
return G1, G2, G3
[docs]def ps(
theta: npt.ArrayLike,
vsvp: Union[float, NDArray],
n: int = 1,
) -> Tuple[NDArray, NDArray, NDArray]:
r"""PS reflection coefficient
Computes the coefficients for the PS approximation
for a set of angles and a constant or variable VS/VP ratio.
Parameters
----------
theta : :obj:`np.ndarray`
Incident angles in degrees
vsvp : :obj:`np.ndarray` or :obj:`float`
:math:`V_S/V_P` ratio
n : :obj:`int`, optional
Number of samples (if ``vsvp`` is a scalar)
Returns
-------
G1 : :obj:`np.ndarray`
First coefficient for VP :math:`[n_{\theta} \times n_\text{vsvp}]`.
Since the PS reflection at zero angle is zero, this value is not used and is
only available to ensure function signature compatibility with other
linearization routines.
G2 : :obj:`np.ndarray`
Second coefficient for VS :math:`[n_{\theta} \times n_\text{vsvp}]`
G3 : :obj:`np.ndarray`
Third coefficient for density :math:`[n_{\theta} \times n_\text{vsvp}]`
Notes
-----
The approximation in [1]_ is used to compute the PS
reflection coefficient as linear combination of contrasts in
:math:`V_P`, :math:`V_S`, and :math:`\rho.` More specifically:
.. math::
R(\theta) = G_2(\theta) \frac{\Delta V_S}{\bar{V_S}} + G_3(\theta)
\frac{\Delta \rho}{\overline{\rho}}
where
.. math::
\begin{align}
G_2(\theta) &= \tan \frac{\theta}{2} \left\{4 (V_S/V_P)^2 \sin^2 \theta
- 4(V_S/V_P) \cos \theta \cos \phi \right\},\\
G_3(\theta) &= -\tan \frac{\theta}{2} \left\{1 - 2 (V_S/V_P)^2 \sin^2 \theta +
2(V_S/V_P) \cos \theta \cos \phi\right\},\\
\frac{\Delta V_S}{\overline{V_S}} &= 2 \frac{V_{S,2}-V_{S,1}}{V_{S,2}+V_{S,1}},\\
\frac{\Delta \rho}{\overline{\rho}} &= 2 \frac{\rho_2-\rho_1}{\rho_2+\rho_1}.
\end{align}
Note that :math:`\theta` is the P-incidence angle whilst :math:`\phi` is
the S-reflected angle which is computed using Snell's law and the average
:math:`V_S/V_P` ratio.
.. [1] Xu, Y., and Bancroft, J.C., "Joint AVO analysis of PP and PS
seismic data", CREWES Report, vol. 9. 1997.
"""
ncp = get_array_module(theta)
theta = ncp.deg2rad(theta)
vsvp = vsvp * np.ones(n) if not isinstance(vsvp, np.ndarray) else vsvp
theta = theta[:, np.newaxis] if vsvp.size > 1 else theta
vsvp = vsvp[:, np.newaxis].T if vsvp.size > 1 else vsvp
phi = np.arcsin(vsvp * np.sin(theta))
# G1 = 0.0 * np.sin(theta) + 0 * vsvp
# G2 = (np.tan(phi) / vsvp) * (4 * np.sin(phi) ** 2 - 4 * vsvp * np.cos(theta) * np.cos(phi)) + 0 * vsvp
# G3 = -((np.tan(phi)) / (2 * vsvp)) * (1 + 2 * np.sin(phi) - 2 * vsvp * np.cos(theta) * np.cos(phi)) + 0 * vsvp
G1 = 0.0 * np.sin(theta) + 0 * vsvp
G2 = (np.tan(phi) / 2) * (
4 * (vsvp * np.sin(phi)) ** 2 - 4 * vsvp * np.cos(theta) * np.cos(phi)
) + 0 * vsvp
G3 = (
-(np.tan(phi) / 2)
* (1 - 2 * (vsvp * np.sin(phi)) ** 2 + 2 * vsvp * np.cos(theta) * np.cos(phi))
+ 0 * vsvp
)
return G1, G2, G3
[docs]class AVOLinearModelling(LinearOperator):
r"""AVO Linearized modelling.
Create operator to be applied to a combination of elastic parameters
for generation of seismic pre-stack reflectivity.
Parameters
----------
theta : :obj:`np.ndarray`
Incident angles in degrees
vsvp : :obj:`np.ndarray` or :obj:`float`
:math:`V_S/V_P` ratio
nt0 : :obj:`int`, optional
Number of samples (if ``vsvp`` is a scalar)
spatdims : :obj:`int` or :obj:`tuple`, optional
Number of samples along spatial axis (or axes)
(``None`` if only one dimension is available)
linearization : `{"akirich", "fatti", "PS"}`, optional
* "akirich": Aki-Richards. See :py:func:`pylops.avo.avo.akirichards`.
* "fatti": Fatti. See :py:func:`pylops.avo.avo.fatti`.
* "PS": PS. See :py:func:`pylops.avo.avo.ps`.
dtype : :obj:`str`, optional
Type of elements in input array.
name : :obj:`str`, optional
.. versionadded:: 2.0.0
Name of operator (to be used by :func:`pylops.utils.describe.describe`)
Attributes
----------
shape : :obj:`tuple`
Operator shape
explicit : :obj:`bool`
Operator contains a matrix that can be solved explicitly
(``True``) or not (``False``)
Raises
------
NotImplementedError
If ``linearization`` is not an implemented linearization
Notes
-----
The AVO linearized operator performs a linear combination of three
(or two) elastic parameters arranged in input vector :math:`\mathbf{m}`
of size :math:`n_{t_0} \times N` to create the so-called seismic
reflectivity:
.. math::
r(t, \theta, x, y) = \sum_{i=1}^N G_i(t, \theta) m_i(t, x, y) \qquad
\forall \,t,\theta
where :math:`N=2,\, 3`. Note that the reflectivity can be in 1d, 2d or 3d
and ``spatdims`` contains the dimensions of the spatial axis (or axes)
:math:`x` and :math:`y`.
"""
def __init__(
self,
theta: NDArray,
vsvp: Union[float, NDArray] = 0.5,
nt0: int = 1,
spatdims: Optional[Union[int, Tuple[int]]] = None,
linearization: str = "akirich",
dtype: DTypeLike = "float64",
name: str = "A",
) -> None:
self.ncp = get_array_module(theta)
self.nt0 = nt0 if not isinstance(vsvp, self.ncp.ndarray) else len(vsvp)
self.ntheta = len(theta)
self.spatdims = () if spatdims is None else _value_or_sized_to_tuple(spatdims)
# Compute AVO coefficients
if linearization == "akirich":
Gs = akirichards(theta, vsvp, n=self.nt0)
elif linearization == "fatti":
Gs = fatti(theta, vsvp, n=self.nt0)
elif linearization == "ps":
Gs = ps(theta, vsvp, n=self.nt0)
else:
logging.error("%s not an available " "linearization...", linearization)
raise NotImplementedError(
"%s not an available linearization..." % linearization
)
self.npars = len(Gs)
dims: Tuple[int, ...] = (self.nt0, self.npars)
dimsd: Tuple[int, ...] = (self.nt0, self.ntheta)
if spatdims is not None:
dims += self.spatdims
dimsd += self.spatdims
super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, name=name)
self.G = self.ncp.concatenate([gs.T[:, self.ncp.newaxis] for gs in Gs], axis=1)
# add dimensions to G to account for horizonal axes
for _ in range(len(self.spatdims)):
self.G = self.G[..., np.newaxis]
@reshaped
def _matvec(self, x: NDArray) -> NDArray:
return self.ncp.sum(self.G * x[:, :, self.ncp.newaxis], axis=1)
@reshaped
def _rmatvec(self, x: NDArray) -> NDArray:
return self.ncp.sum(self.G * x[:, self.ncp.newaxis], axis=2)
