Basis sets¶
- class mumott.methods.basis_sets.SphericalHarmonics(probed_coordinates=None, ell_max=0, enforce_friedel_symmetry=True, **kwargs)[source]¶
Basis set class for spherical harmonics, the canonical representation of polynomials on the unit sphere and a simple way of representing band-limited spherical functions which allows for easy computations of statistics and is suitable for analyzing certain symmetries.
- Parameters
ell_max (int) – The bandlimit of the spherical functions that you want to be able to represent. A good rule of thumb is that
ell_max
should not exceed the number of detector segments minus 1.probed_coordinates (ProbedCoordinates) – Optional. A container with the coordinates on the sphere probed at each detector segment by the experimental method. Its construction from the system geometry is method-dependent. By default, an empty instance of
ProbedCoordinates
is created.enforce_friedel_symmetry (bool) – If set to
True
, Friedel symmetry will be enforced, using the assumption that points on opposite sides of the sphere are equivalent. This results in only evenell
being used.kwargs –
Miscellaneous arguments which relate to segment integrations can be passed as keyword arguments:
- integration_mode
Mode to integrate line segments on the reciprocal space sphere. Possible options are
'simpson'
,'midpoint'
,'romberg'
,'trapezoid'
.'simpson'
,'trapezoid'
, and'romberg'
use adaptive integration with the respective quadrature rule fromscipy.integrate
.'midpoint'
uses a single mid-point approximation of the integral. Default value is'simpson'
.- n_integration_starting_points
Number of points used in the first iteration of the adaptive integration. The number increases by the rule
N
←2 * N - 1
for each iteration. Default value is 3.- integration_tolerance
Tolerance for the maximum relative error between iterations before the integral is considered converged. Default is
1e-5
.- integration_maxiter
Maximum number of iterations. Default is
10
.
- property csr_representation: tuple¶
The projection matrix as a stack of sparse matrices in CSR representation as a tuple. The information in the tuple consists of the 3 dense matrices making up the representation, in the order
(pointers, indices, data)
.
- property ell_indices: ndarray[Any, dtype[int]]¶
The
ell
associated with each coefficient and its corresponding spherical harmonic. Updated whenell_max
changes.Notes
The word
ell
is used to represent the cursive small L, also written \(\ell\), often used as an index for the degree of the Legendre polynomial in the definition of the spherical harmonics.
- property ell_max: int¶
The maximum
ell
used to represent spherical functions.Notes
The word
ell
is used to represent the cursive small L, also written \(\ell\), often used as an index for the degree of the Legendre polynomial in the definition of the spherical harmonics.
- property emm_indices: ndarray[Any, dtype[int]]¶
The
emm
associated with each coefficient and its corresponding spherical harmonic. Updated whenell_max
changes.Notes
For consistency with
ell_indices
, and to avoid visual confusion with other letters,emm
is used to represent the index commonly written \(m\) in mathematical notation, the frequency of the sine-cosine parts of the spherical harmonics, often called the spherical harmonic order.
- forward(coefficients, indices=None)[source]¶
Carries out a forward computation of projections from spherical harmonic space into detector space, for one or several tomographic projections.
- Parameters
coefficients (
array
) – An array of coefficients, of arbitrary shape so long as the last axis has the same size asell_indices
, and ifindices
is None or greater than one, the first axis should have the same length asindices
indices (
Optional
[array
]) – Optional. Indices of the tomographic projections for which the forward computation is to be performed. IfNone
, the forward computation will be performed for all projections.
- Return type
array
- Returns
An array of values on the detector corresponding to the
coefficients
given. Ifindices
contains exactly one index, the shape is(coefficients.shape[:-1], J)
whereJ
is the number of detector segments. Ifindices
isNone
or contains several indices, the shape is(N, coefficients.shape[1:-1], J)
whereN
is the number of tomographic projections for which the computation is performed.
Notes
The assumption is made in this implementation that computations over several indices act on sets of images from different projections. For special usage where multiple projections of entire fields is desired, it may be better to use
projection_matrix
directly. This also applies togradient()
.
- generate_map(coefficients, resolution_in_degrees=5, map_half_sphere=True)¶
Generate a (theta, phi) map of the function modeled by the input coefficients. If
map_half_sphere=True
(default) a map of only the z>0 half sphere is returned.- Parameters
coefficients (
ndarray
[Any
,dtype
[float
]]) – One dimensional numpy array with lengthlen(self)
containing the coefficients of the function to be plotted.resolution_in_degrees (
int
) – The resoution of the map in degrees. The map uses eqidistant lines in longitude and latitude.map_half_sphere (
bool
) – If True returns a map of the z>0 half sphere.
- Return type
tuple
[ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]]- Returns
map_intensity – Intensity values of the map.
map_theta – Polar cooridnates of the map.
map_phi – Azimuthal coordinates of the map.
- get_covariances(u, v, resolve_spectrum=False)[source]¶
Returns the covariances of the spherical functions represented by two coefficient arrays.
- Parameters
u (
array
) – The first coefficient array, of arbitrary shape except its last dimension must be the same length as the length ofell_indices
.v (
array
) – The second coefficient array, of the same shape asu
.resolve_spectrum (
bool
) – Optional. Whether to resolve the product according to each frequency band, given by the coefficients of eachell
inell_indices
. Default value isFalse
.
- Return type
array
- Returns
An array of the covariances of the spherical functions represented by
u
andv
. Has the shape(u.shape[:-1])
if resolve_spectrum isFalse
, and(u.shape[:-1] + (ell_max // 2 + 1,))
if resolve_spectrum isTrue
, whereell_max
isSphericalHarmonics.ell_max
.
Notes
Calling this function is equivalent to calling
get_inner_product()
withspectral_moments=np.unique(ell_indices[ell_indices > 0])
whereell_indices
isSphericalHarmonics.ell_indices
. See the note toget_inner_product()
for mathematical details.
- get_inner_product(u, v, resolve_spectrum=False, spectral_moments=None)[source]¶
Retrieves the inner product of two coefficient arrays.
Notes
The canonical inner product in a spherical harmonic representation is \(\sum_\ell N(\ell) \sum_m u_m^\ell v_m^\ell\), where \(N(\ell)\) is a normalization constant (which is unity for the \(4\pi\) normalization). This inner product is a rotational invariant. The rotational invariance also holds for any partial sums over \(\ell\). One can define a function of \(\ell\) that returns such products, namely \(S(\ell, u, v) = N(\ell)\sum_m u_m^\ell v_m^\ell\), called the spectral power function. The sum \(\sum_{\ell = 1}S(\ell)\) is equal to the covariance of the band-limited spherical functions represented by \(u\) and \(v\), and each \(S(\ell, u, v)\) is the contribution to the covariance of the band \(\ell\). See also the SHTOOLS documentation for an excellent overview of this.
- Parameters
u (
array
) – The first coefficient array, of arbitrary shape and dimension, except the last dimension must be the same as the length ofell_indices
.v (
array
) – The second coefficient array, of the same shape asu
.resolve_spectrum (
bool
) – Optional. Whether to resolve the product according to each frequency band, given by the coefficients of eachell
inell_indices
. Defaults toFalse
, which means that the sum of every component of the spectrum is returned. IfTrue
, components are returned in order of ascendingell
. Theell
included in the spectrum depends onspectral_moments
.spectral_moments (
Optional
[List
[int
]]) – Optional. List of particular values ofell
to calculate the inner product for. Defaults toNone
, which is identical to including all values ofell
in the calculation. Ifspectral_moments
contains all nonzero values ofell
andresolve_spectrum
isFalse
, the covariance ofv
andu
will be calculated (the sum of the inner product over all non-zeroell
Ifresolve_spectrum
isTrue
, the covariance per ell inspectral_moments
, will be calculated, i.e., the inner products will not be summed over.
- Return type
array
- Returns
An array of the inner products of the spherical functions represented by
u
andv
. Has the shape(u.shape[:-1])
ifresolve_spectrum
isFalse
,(u.shape[:-1] + (ell_max // 2 + 1,))
ifresolve_spectrum
isTrue
andspectral_moments
isNone
, and finally the shape(u.shape[:-1] + (np.unique(spectral_moments).size,))
ifresolve_spectrum
isTrue
andspectral_moments
is a list of integers found inell_indices
- get_output(coefficients)[source]¶
Returns a
ReconstructionDerivedQuantities
instance of output data for a given array of basis set coefficients.- Parameters
coefficients (
array
) – An array of coefficients of arbitrary shape and dimensions, except its last dimension must be the same length as thelen
of this instance. Computations only operate over the last axis ofcoefficients
, so derived properties in the output will have the shape(*coefficients.shape[:-1], ...)
.- Return type
- Returns
ReconstructionDerivedQuantities
containing a number of quantities that have been computed from the spherical functions represented by the input coefficients.
- get_spherical_harmonic_coefficients(coefficients, ell_max=None)[source]¶
Convert a set of spherical harmonics coefficients to a different
ell_max
by either zero-padding or truncation and return the result.
- gradient(coefficients, indices=None)[source]¶
Carries out a gradient computation of projections from spherical harmonic space into detector space, for one or several tomographic projections.
- Parameters
coefficients (
array
) – An array of coefficients (or residuals) of arbitrary shape so long as the last axis has the same size as the number of detector segments.indices (
Optional
[array
]) – Optional. Indices of the tomographic projections for which the gradient computation is to be performed. IfNone
, the gradient computation will be performed for all projections.
- Return type
array
- Returns
An array of gradient values based on the coefficients given. If
indices
contains exactly one index, the shape is(coefficients.shape[:-1], J)
whereJ
is the number of detector segments. If indices isNone
or contains several indices, the shape is(N, coefficients.shape[1:-1], J)
whereN
is the number of tomographic projections for which the computation is performed.
Notes
When solving an inverse problem, one should not to attempt to optimize the coefficients directly using the
gradient
one obtains by applying this method to the data. Instead, one must either take the gradient of the residual between theforward()
computation of the coefficients and the data. Alternatively one can apply both the forward and the gradient computation to the coefficients to be optimized, and the gradient computation to the data, and treat the residual of the two as the gradient of the optimization coefficients. The approaches are algebraically equivalent, but one may be more efficient than the other in some circumstances.
- property integration_mode: str¶
Mode of integration for calculating projection matrix. Accepted values are
'simpson'
,'romberg'
,'trapezoid'
, and'midpoint'
.
- property projection_matrix: array¶
The matrix used to project spherical functions from the unit sphere onto the detector. If
v
is a vector of spherical harmonic coefficients, andM
is theprojection_matrix
, thenM @ v
gives the corresponding values on the detector segments associated with each projection.M[i] @ v
gives the values on the detector segments associated with projectioni
.If
r
is a residual between a projection from spherical to detector space and data from projectioni
, thenM[i].T @ r
gives the associated gradient in spherical harmonic space.
- class mumott.methods.basis_sets.TrivialBasis(channels=1)[source]¶
Basis set class for the trivial basis, i.e., the identity basis. This can be used as a scaffolding class when implementing, e.g., scalar tomography, as it implements all the necessary functionality to qualify as a
BasisSet
.- Parameters
channels (
int
) – Number of channels in the last index. Default is1
. For scalar data, the default value of1
is appropriate. For any other use-case, where the representation on the sphere and the representation in detector space are equivalent, such as reconstructing scalars of multiple q-ranges at once, a different number of channels can be set.
- property csr_representation: tuple¶
The projection matrix as a stack of sparse matrices in CSR representation as a tuple. The information in the tuple consists of the 3 dense matrices making up the representation, in the order
(pointers, indices, data)
.
- forward(coefficients, *args, **kwargs)[source]¶
Returns the provided coefficients with no modification.
- Parameters
coefficients (
ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]) – An array of coefficients, of arbitrary shape, except the last index must specify the same number of channels as was specified for this basis.- Return type
ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]- Returns
The provided :attr`coefficients` with no modification.
Notes
The
args
andkwargs
are ignored, but included for compatibility with methods that input other arguments.
- generate_map(coefficients, resolution_in_degrees=5, map_half_sphere=True)¶
Generate a (theta, phi) map of the function modeled by the input coefficients. If
map_half_sphere=True
(default) a map of only the z>0 half sphere is returned.- Parameters
coefficients (
ndarray
[Any
,dtype
[float
]]) – One dimensional numpy array with lengthlen(self)
containing the coefficients of the function to be plotted.resolution_in_degrees (
int
) – The resoution of the map in degrees. The map uses eqidistant lines in longitude and latitude.map_half_sphere (
bool
) – If True returns a map of the z>0 half sphere.
- Return type
tuple
[ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]]- Returns
map_intensity – Intensity values of the map.
map_theta – Polar cooridnates of the map.
map_phi – Azimuthal coordinates of the map.
- get_inner_product(u, v)[source]¶
Retrieves the inner product of two coefficient arrays, that is to say, the sum-product over the last axis.
- get_output(coefficients)[source]¶
Returns a dictionary of output data for a given array of coefficients.
- Parameters
coefficients (
ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]) – An array of coefficients of arbitrary shape and dimension. Computations only operate over the last axis ofcoefficents
, so derived properties in the output will have the shape(*coefficients.shape[:-1], ...)
.- Return type
- Returns
A dictionary containing a dictionary with the field
basis_set
.
Notes
In detail, the dictionary under the key
basis_set
contains:- basis_set
- name
The name of the basis set, i.e.,
'TrivialBasis'
- coefficients
A copy of
coefficients
.- projection_matrix
The identity matrix of the same size as the number of chanenls.
- get_spherical_harmonic_coefficients(coefficients, ell_max=None)[source]¶
Convert a set of spherical harmonics coefficients to a different
ell_max
by either zero-padding or truncation and return the result.
- gradient(coefficients, *args, **kwargs)[source]¶
Returns the provided coefficients with no modification.
- Parameters
coefficients (
ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]) – An array of coefficients of arbitrary shape except the last index must specify the same number of channels as was specified for this basis.- Return type
ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]- Returns
The provided :attr`coefficients` with no modification.
Notes
The
args
andkwargs
are ignored, but included for compatibility with methods that input other argumetns.
- property integration_mode: str¶
Mode of integration for calculating projection matrix. Accepted values are
'simpson'
,'romberg'
,'trapezoid'
, and'midpoint'
.
- property probed_coordinates: ProbedCoordinates¶
- property projection_matrix¶
The identity matrix of the same rank as the number of channels specified.
- class mumott.methods.basis_sets.GaussianKernels(probed_coordinates=None, grid_scale=4, kernel_scale_parameter=1.0, enforce_friedel_symmetry=True, **kwargs)[source]¶
Basis set class for gaussian kernels, a simple local representation on the sphere. The kernels follow a pseudo-even distribution similar to that described by Y. Kurihara in 1965, except with offsets added at the poles.
Notes
The Gaussian kernel at location \(\rho_i\) is given by
\[N_i \exp\left[ -\frac{1}{2} \left(\frac{d(\rho_i, r)}{\sigma}\right)^2 \right]\]\[\sigma = \frac{\nu \pi}{2 (g + 1)}\]where \(\nu\) is the kernel scale parameter and \(g\) is the grid scale, and
\[d(\rho, r) = \arctan_2(\Vert \rho \times r \Vert, \rho \cdot r),\]that is, the great circle distance from the kernel location \(\rho\) to the probed location \(r\). If Friedel symmetry is assumed, the expression is instead
\[d(\rho, r) = \arctan_2(\Vert \rho \times r \Vert, \vert \rho \cdot r \vert)\]The normalization factor \(\rho_i\) is given by
\[N_i = \sum_j \exp\left[ -\frac{1}{2} \left( \frac{d(\rho_i, \rho_j)}{\sigma} \right)^2 \right]\]where the sum goes over the coordinates of all grid points. This leads to an approximately even spherical function, such that a set of coefficients which are all equal is approximately isotropic, to the extent possible with respect to restrictions imposed by grid resolution and scale parameter.
- Parameters
probed_coordinates (ProbedCoordinates) – Optional. A container with the coordinates on the sphere probed at each detector segment by the experimental method. Its construction from the system geometry is method-dependent. By default, an empty instance of
mumott.ProbedCoordinates
is created.grid_scale (int) – The size of the coordinate grid on the sphere. Denotes the number of azimuthal rings between the pole and the equator, where each ring has between
2
and2 * grid_scale
points along the azimuth.kernel_scale_parameter (float) – The scale parameter of the kernel in units of \(\frac{\pi}{2 (g + 1)}\), where \(g\) is
grid_scale
.enforce_friedel_symmetry (bool) – If set to
True
, Friedel symmetry will be enforced, using the assumption that points on opposite sides of the sphere are equivalent.kwargs –
Miscellaneous arguments which relate to segment integrations can be passed as keyword arguments:
- integration_mode
Mode to integrate line segments on the reciprocal space sphere. Possible options are
'simpson'
,'midpoint'
,'romberg'
,'trapezoid'
.'simpson'
,'trapezoid'
, and'romberg'
use adaptive integration with the respective quadrature rule fromscipy.integrate
.'midpoint'
uses a single mid-point approximation of the integral. Default value is'simpson'
.- n_integration_starting_points
Number of points used in the first iteration of the adaptive integration. The number increases by the rule
N
←2 * N - 1
for each iteration. Default value is 3.- integration_tolerance
Tolerance for the maximum relative error between iterations before the integral is considered converged. Default is
1e-5
.- integration_maxiter
Maximum number of iterations. Default is
10
.- enforce_sparsity
If
True
, makes matrix sparse by limiting the number of basis set elements that can map to each segment. Default isFalse
.- sparsity_count
Number of basis set elements that can map to each segment, if
enforce_sparsity
is set toTrue
. Default is3
.
- property csr_representation: tuple¶
The projection matrix as a stack of sparse matrices in CSR representation as a tuple. The information in the tuple consists of the 3 dense matrices making up the representation, in the order
(pointers, indices, data)
.
- property enforce_friedel_symmetry: bool¶
If
True
, Friedel symmetry is enforced, i.e., the point \(-r\) is treated as equivalent to \(r\).
- forward(coefficients, indices=None)[source]¶
Carries out a forward computation of projections from Gaussian kernel space into detector space, for one or several tomographic projections.
- Parameters
coefficients (
ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]) – An array of coefficients, of arbitrary shape so long as the last axis has the same size askernel_scale_parameter
, and ifindices
isNone
or greater than one, the first axis should have the same length asindices
indices (
Optional
[ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]]) – Optional. Indices of the tomographic projections for which the forward computation is to be performed. IfNone
, the forward computation will be performed for all projections.
- Return type
ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]- Returns
An array of values on the detector corresponding to the
coefficients
given. Ifindices
contains exactly one index, the shape is(coefficients.shape[:-1], J)
whereJ
is the number of detector segments. Ifindices
isNone
or contains several indices, the shape is(N, coefficients.shape[1:-1], J)
whereN
is the number of tomographic projections for which the computation is performed.
Notes
The assumption is made in this implementation that computations over several indices act on sets of images from different projections. For special usage where multiple projections of entire fields are desired, it may be better to use
projection_matrix
directly. This also applies togradient()
.
- generate_map(coefficients, resolution_in_degrees=5, map_half_sphere=True)¶
Generate a (theta, phi) map of the function modeled by the input coefficients. If
map_half_sphere=True
(default) a map of only the z>0 half sphere is returned.- Parameters
coefficients (
ndarray
[Any
,dtype
[float
]]) – One dimensional numpy array with lengthlen(self)
containing the coefficients of the function to be plotted.resolution_in_degrees (
int
) – The resoution of the map in degrees. The map uses eqidistant lines in longitude and latitude.map_half_sphere (
bool
) – If True returns a map of the z>0 half sphere.
- Return type
tuple
[ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]]- Returns
map_intensity – Intensity values of the map.
map_theta – Polar cooridnates of the map.
map_phi – Azimuthal coordinates of the map.
- get_amplitudes(coefficients, probed_coordinates=None)[source]¶
Computes the amplitudes of the spherical function represented by the provided
coefficients
at theprobed_coordinates
.- Parameters
coefficients (
ndarray
[Any
,dtype
[float
]]) – An array of coefficients of arbitrary shape, provided that the last dimension contains the coefficients for one spherical function.probed_coordinates (
Optional
[ProbedCoordinates
]) – An instance ofmumott.core.ProbedCoordinates
with itsvector
attribute indicating the points of the sphere for which to evaluate the amplitudes.
- Return type
- get_inner_product(u, v)[source]¶
Retrieves the inner product of two coefficient arrays, that is to say, the sum-product over the last axis.
- Parameters
u (
ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]) – The first coefficient array, of arbitrary shape and dimension, so long as the number of coefficients equals the length of thisGaussianKernels
instance.v (
ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]) – The second coefficient array, of the same shape asu
.
- Return type
ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]
- get_output(coefficients)[source]¶
Returns a
ReconstructionDerivedQuantities
instance of output data for a given array of basis set coefficients.- Parameters
coefficients (
ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]) – An array of coefficients of arbitrary shape and dimensions, except its last dimension must be the same length as thelen
of this instance. Computations only operate over the last axis ofcoefficients
, so derived properties in the output will have the shape(*coefficients.shape[:-1], ...)
.- Return type
- Returns
ReconstructionDerivedQuantities
containing a number of quantities that have been computed from the spherical functions represented by the input coefficients.
- get_second_moments(coefficients)[source]¶
Calculate the second moments of the functions described by
coefficients
.- Parameters
coefficients (
ndarray
[Any
,dtype
[float
]]) – An array of coefficients (or residuals) of arbitrary shape so long as the last axis has the same size as the number of detector channels.- Return type
- Returns
Array containing the second moments of the functions described by coefficients, formatted as rank-two tensors with tensor indices in the last 2 dimensions.
- get_spherical_harmonic_coefficients(coefficients, ell_max=None)[source]¶
Computes the spherical harmonic coefficients of the spherical function represented by the provided
coefficients
using a Driscoll-Healy grid.For details on the Driscoll-Healy grid, see the SHTools page for a comprehensive overview.
- Parameters
coefficients (
ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]) – An array of coefficients of arbitrary shape, provided that the last dimension contains the coefficients for one function.ell_max (
Optional
[int
]) – The bandlimit of the spherical harmonic expansion. By default, it is2 * grid_scale
.
- gradient(coefficients, indices=None)[source]¶
Carries out a gradient computation of projections from Gaussian kernel space into detector space for one or several tomographic projections.
- Parameters
coefficients (
ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]) – An array of coefficients (or residuals) of arbitrary shape so long as the last axis has the same size as the number of detector segments.indices (
Optional
[ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]]) – Optional. Indices of the tomographic projections for which the gradient computation is to be performed. IfNone
, the gradient computation will be performed for all projections.
- Return type
ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]- Returns
An array of gradient values based on the
coefficients
given. Ifindices
contains exactly one index, the shape is(coefficients.shape[:-1], J)
whereJ
is the number of detector segments. If indices isNone
or contains several indices, the shape is(N, coefficients.shape[1:-1], J)
whereN
is the number of tomographic projections for which the computation is performed.
Notes
When solving an inverse problem, one should not attempt to optimize the coefficients directly using the gradient one obtains by applying this method to the data. Instead, one must either take the gradient of the residual between the
forward()
computation of the coefficients and the data. Alternatively one can apply both the forward and the gradient computation to the coefficients to be optimized, and the gradient computation to the data, and treat the residual of the two as the gradient of the optimization coefficients. The approaches are algebraically equivalent, but one may be more efficient than the other in some circumstances. However, normally, the projection between detector andGaussianKernel
space is only a small part of the overall computation, so there is typically not much to be gained from optimizing it.
- property grid: Tuple[ndarray[Any, dtype[float]], ndarray[Any, dtype[float]]]¶
Returns the polar and azimuthal angles of the grid used by the basis.
- Returns
A
Tuple
with contents(polar_angle, azimuthal_angle)
, where thepolar angle is defined as \(\arccos(z)\).
- property grid_scale: int¶
The number of azimuthal rings from each pole to the equator in the spherical grid.
- property integration_mode: str¶
Mode of integration for calculating projection matrix. Accepted values are
'simpson'
,'romberg'
,'trapezoid'
, and'midpoint'
.
- property probed_coordinates: ProbedCoordinates¶
- property projection_matrix: ndarray[Any, dtype[_ScalarType_co]]¶
The matrix used to project spherical functions from the unit sphere onto the detector. If
v
is a vector of gaussian kernel coefficients, andM
is theprojection_matrix
, thenM @ v
gives the corresponding values on the detector segments associated with each projection.M[i] @ v
gives the values on the detector segments associated with projectioni
.If
r
is a residual between a projection from Gaussian kernel to detector space and data from projectioni
, thenM[i].T @ r
gives the associated gradient in Gaussian kernel space.
- property projection_matrix_hash: str¶
Returns a hash of
projection_matrix
.
- class mumott.methods.basis_sets.NearestNeighbor(directions, probed_coordinates=None, enforce_friedel_symmetry=True, **kwargs)[source]¶
Basis set class for nearest-neighbor interpolation. Used to construct methods similar to that presented in Schaff et al. (2015). By default this representation is sparse and maps only a single direction on the sphere to each detector segment. This can be changed; see
kwargs
.- Parameters
directions (NDArray[float]) – Two-dimensional Array containing the
N
sensitivity directions with shape(N, 3)
.probed_coordinates (ProbedCoordinates) – Optional. Coordinates on the sphere probed at each detector segment by the experimental method. Its construction from the system geometry is method-dependent. By default, an empty instance of
mumott.ProbedCoordinates
is created.enforce_friedel_symmetry (bool) – If set to
True
, Friedel symmetry will be enforced, using the assumption that points on opposite sides of the sphere are equivalent.kwargs –
Miscellaneous arguments which relate to segment integrations can be passed as keyword arguments:
- integration_mode
Mode to integrate line segments on the reciprocal space sphere. Possible options are
'simpson'
,'midpoint'
,'romberg'
,'trapezoid'
.'simpson'
,'trapezoid'
, and'romberg'
use adaptive integration with the respective quadrature rule fromscipy.integrate
.'midpoint'
uses a single mid-point approximation of the integral. Default value is'simpson'
.- n_integration_starting_points
Number of points used in the first iteration of the adaptive integration. The number increases by the rule
N
←2 * N - 1
for each iteration. Default value is 3.- integration_tolerance
Tolerance for the maximum relative error between iterations before the integral is considered converged. Default is
1e-3
.- integration_maxiter
Maximum number of iterations. Default is
10
.- enforce_sparsity
If
True
, limites the number of basis set elements that can map to each detector segemnt. Default isFalse
.- sparsity_count
If
enforce_sparsity
is set toTrue
, the number of basis set elements that can map to each detector segment. Default value is1
.
- property csr_representation: tuple¶
The projection matrix as a stack of sparse matrices in CSR representation as a tuple. The information in the tuple consists of the 3 dense matrices making up the representation, in the order
(pointers, indices, data)
.
- property enforce_friedel_symmetry: bool¶
If
True
, Friedel symmetry is enforced, i.e., the point \(-r\) is treated as equivalent to \(r\).
- find_nearest_neighbor_index(probed_directions)[source]¶
Caluculate the nearest neighbor sensitivity directions for an array of x-y-z vectors.
- forward(coefficients, indices=None)[source]¶
Carries out a forward computation of projections from reciprocal space modes to detector channels, for one or several tomographic projections.
- Parameters
coefficients (
ndarray
[Any
,dtype
[float
]]) – An array of coefficients, of arbitrary shape so long as the last axis has the same size as this basis set.indices (
Optional
[ndarray
[Any
,dtype
[int
]]]) – Optional. Indices of the tomographic projections for which the forward computation is to be performed. IfNone
, the forward computation will be performed for all projections.
- Return type
- Returns
An array of values on the detector corresponding to the
coefficients
given. Ifindices
contains exactly one index, the shape is(coefficients.shape[:-1], J)
whereJ
is the number of detector segments. Ifindices
isNone
or contains several indices, the shape is(N, coefficients.shape[1:-1], J)
whereN
is the number of tomographic projections for which the computation is performed.
- generate_map(coefficients, resolution_in_degrees=5, map_half_sphere=True)¶
Generate a (theta, phi) map of the function modeled by the input coefficients. If
map_half_sphere=True
(default) a map of only the z>0 half sphere is returned.- Parameters
coefficients (
ndarray
[Any
,dtype
[float
]]) – One dimensional numpy array with lengthlen(self)
containing the coefficients of the function to be plotted.resolution_in_degrees (
int
) – The resoution of the map in degrees. The map uses eqidistant lines in longitude and latitude.map_half_sphere (
bool
) – If True returns a map of the z>0 half sphere.
- Return type
tuple
[ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]]- Returns
map_intensity – Intensity values of the map.
map_theta – Polar cooridnates of the map.
map_phi – Azimuthal coordinates of the map.
- get_amplitudes(coefficients, probed_directions)[source]¶
Calculate function values of an array of coefficients.
- Parameters
- Return type
- Returns
Array with function values. The shape of the array is
(*coefficients.shape[:-1], *probed_directions.shape[:-1])
.
- get_function_values(probed_directions)[source]¶
Calculate the value of the basis functions from an array of x-y-z vectors.
- Parameters
probed_directions (
ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]) – Array with length 3 along its last axis- Return type
- Returns
Array with same shape as input array except for the last axis, which now has length
N
, i.e., the number of sensitivity directions.
- get_output(coefficients)[source]¶
Returns a
ReconstructionDerivedQuantities
instance of output data for a given array of basis set coefficients.- Parameters
coefficients (
ndarray
[Any
,dtype
[TypeVar
(_ScalarType_co
, bound=generic
, covariant=True)]]) – An array of coefficients of arbitrary shape and dimensions, except its last dimension must be the same length as thelen
of this instance. Computations only operate over the last axis ofcoefficients
, so derived properties in the output will have the shape(*coefficients.shape[:-1], ...)
.- Return type
- Returns
ReconstructionDerivedQuantities
containing a number of quantities that have been computed from the spherical functions represented by the input coefficients.
- get_second_moments(coefficients)[source]¶
Calculate the second moments of the functions described by
coefficients
.- Parameters
coefficients (
ndarray
[Any
,dtype
[float
]]) – An array of coefficients (or residuals) of arbitrary shape so long as the last axis has the same size as the number of detector channels.- Return type
- Returns
Array containing the second moments of the functions described by coefficients, formatted as rank-two tensors with tensor indices in the last 2 dimensions.
- get_spherical_harmonic_coefficients(coefficients, ell_max=None)[source]¶
Computes and rturns the spherical harmonic coefficients of the spherical function represented by the provided
coefficients
using a Driscoll-Healy grid.For details on the Driscoll-Healy grid, see the SHTools page for a comprehensive overview.
- get_sub_geometry(direction_index, geometry, data_container=None)[source]¶
Create and return a geometry object corresponding to a scalar tomography problem for scattering along the sensitivity direction with index
direction_index
. If optionally amumott.DataContainer
is provided, the sinograms and weights for this scalar tomography problem will alse be returned.Used for an implementation of the algorithm descibed in [Schaff2015].
- gradient(coefficients, indices=None)[source]¶
Carries out a gradient computation of projections of projections from reciprocal space modes to detector channels, for one or several tomographic projections.
- Parameters
coefficients (
ndarray
[Any
,dtype
[float
]]) – An array of coefficients (or residuals) of arbitrary shape so long as the last axis has the same size as the number of detector channels.indices (
Optional
[ndarray
[Any
,dtype
[int
]]]) – Optional. Indices of the tomographic projections for which the gradient computation is to be performed. IfNone
, the gradient computation will be performed for all projections.
- Return type
- Returns
An array of gradient values based on the
coefficients
given. Ifindices
contains exactly one index, the shape is(coefficients.shape[:-1], J)
whereJ
is the number of detector segments. If indices isNone
or contains several indices, the shape is(N, coefficients.shape[1:-1], J)
whereN
is the number of tomographic projections for which the computation is performed.
- property grid: Tuple[ndarray[Any, dtype[float]], ndarray[Any, dtype[float]]]¶
Returns the polar and azimuthal angles of the grid used by the basis.
- Returns
A
Tuple
with contents(polar_angle, azimuthal_angle)
, where thepolar angle is defined as \(\arccos(z)\).
- property integration_mode: str¶
Mode of integration for calculating projection matrix. Accepted values are
'simpson'
,'romberg'
,'trapezoid'
, and'midpoint'
.
- property probed_coordinates: ProbedCoordinates¶
- property projection_matrix: ndarray[Any, dtype[_ScalarType_co]]¶
The matrix used to project spherical functions from the unit sphere onto the detector. If
v
is a vector of gaussian kernel coefficients, andM
is theprojection_matrix
, thenM @ v
gives the corresponding values on the detector segments associated with each projection.M[i] @ v
gives the values on the detector segments associated with projectioni
.
- property projection_matrix_hash: str¶
Returns a hash of
projection_matrix
.