Coverage for local_installation_linux/mumott/methods/basis_sets/gaussian_kernels.py: 90%
233 statements
« prev ^ index » next coverage.py v7.3.2, created at 2025-05-05 21:21 +0000
1import logging
2from typing import Any, Iterator, Tuple
4import numpy as np
5from numpy.typing import NDArray
7from mumott import ProbedCoordinates, SphericalHarmonicMapper
8from mumott.core.hashing import list_to_hash
9from mumott.methods.utilities.tensor_operations import (framewise_contraction,
10 framewise_contraction_transpose)
11from mumott.output_handling.reconstruction_derived_quantities import\
12 (ReconstructionDerivedQuantities, get_sorted_eigenvectors)
13from .base_basis_set import BasisSet
16logger = logging.getLogger(__name__)
19class GaussianKernels(BasisSet):
20 r""" Basis set class for gaussian kernels, a simple local representation on the sphere.
21 The kernels follow a pseudo-even distribution similar to that described by
22 Y. Kurihara `in 1965 <https://doi.org/10.1175/1520-0493%281965%29093%3C0399%3ANIOTPE%3E2.3.CO%3B2>`_,
23 except with offsets added at the poles.
25 Notes
26 -----
27 The Gaussian kernel at location :math:`\rho_i` is given by
29 .. math ::
30 N_i \exp\left[ -\frac{1}{2} \left(\frac{d(\rho_i, r)}{\sigma}\right)^2 \right]
32 .. math ::
33 \sigma = \frac{\nu \pi}{2 (g + 1)}
35 where :math:`\nu` is the kernel scale parameter and :math:`g` is the grid scale, and
37 .. math ::
38 d(\rho, r) = \arctan_2(\Vert \rho \times r \Vert, \rho \cdot r),
40 that is, the great circle distance from the kernel location :math:`\rho` to the
41 probed location :math:`r`. If Friedel symmetry is assumed, the expression is instead
43 .. math ::
44 d(\rho, r) = \arctan_2(\Vert \rho \times r \Vert, \vert \rho \cdot r \vert)
46 The normalization factor :math:`\rho_i` is given by
48 .. math ::
49 N_i = \sum_j \exp\left[ -\frac{1}{2} \left( \frac{d(\rho_i, \rho_j)}{\sigma} \right)^2 \right]
51 where the sum goes over the coordinates of all grid points. This leads to an
52 approximately even spherical function, such that a set of coefficients which are all equal
53 is approximately isotropic, to the extent possible with respect to restrictions
54 imposed by grid resolution and scale parameter.
56 Parameters
57 ----------
58 probed_coordinates : ProbedCoordinates
59 Optional. A container with the coordinates on the sphere probed at each detector segment by the
60 experimental method. Its construction from the system geometry is method-dependent.
61 By default, an empty instance of :class:`mumott.ProbedCoordinates` is created.
62 grid_scale : int
63 The size of the coordinate grid on the sphere. Denotes the number of azimuthal rings between the
64 pole and the equator, where each ring has between ``2`` and ``2 * grid_scale`` points
65 along the azimuth.
66 kernel_scale_parameter : float
67 The scale parameter of the kernel in units of :math:`\frac{\pi}{2 (g + 1)}`, where
68 :math:`g` is ``grid_scale``.
69 enforce_friedel_symmetry : bool
70 If set to ``True``, Friedel symmetry will be enforced, using the assumption that points
71 on opposite sides of the sphere are equivalent.
72 kwargs
73 Miscellaneous arguments which relate to segment integrations can be
74 passed as keyword arguments:
76 integration_mode
77 Mode to integrate line segments on the reciprocal space sphere. Possible options are
78 ``'simpson'``, ``'midpoint'``, ``'romberg'``, ``'trapezoid'``.
79 ``'simpson'``, ``'trapezoid'``, and ``'romberg'`` use adaptive
80 integration with the respective quadrature rule from ``scipy.integrate``.
81 ``'midpoint'`` uses a single mid-point approximation of the integral.
82 Default value is ``'simpson'``.
83 n_integration_starting_points
84 Number of points used in the first iteration of the adaptive integration.
85 The number increases by the rule ``N`` ← ``2 * N - 1`` for each iteration.
86 Default value is 3.
87 integration_tolerance
88 Tolerance for the maximum relative error between iterations before the integral
89 is considered converged. Default is ``1e-5``.
90 integration_maxiter
91 Maximum number of iterations. Default is ``10``.
92 enforce_sparsity
93 If ``True``, makes matrix sparse by limiting the number of basis set elements
94 that can map to each segment. Default is ``False``.
95 sparsity_count
96 Number of basis set elements that can map to each segment,
97 if ``enforce_sparsity`` is set to ``True``. Default is ``3``.
98 """
99 def __init__(self,
100 probed_coordinates: ProbedCoordinates = None,
101 grid_scale: int = 4,
102 kernel_scale_parameter: float = 1.,
103 enforce_friedel_symmetry: bool = True,
104 **kwargs):
105 super().__init__(probed_coordinates, **kwargs)
106 self._probed_coordinates_hash = hash(self.probed_coordinates)
107 self._grid_scale = grid_scale
108 self._kernel_scale_parameter = kernel_scale_parameter
109 self._enforce_friedel_symmetry = enforce_friedel_symmetry
110 self._projection_matrix = self._get_integrated_projection_matrix()
112 def _get_kurihara_mesh(self, N) -> Tuple[NDArray, NDArray]:
113 phi = []
114 theta = []
115 for i in np.arange(N, -1, -1):
116 for j in np.arange(0, (2 * (i + 0.5))):
117 phi.append((j + 0.5) / (2 * (i + 0.5)) * np.pi)
118 phi.append((j + 0.5) / (2 * (i + 0.5)) * -np.pi)
119 theta.append((i + 0.5) / (N + 1) * np.pi / 2)
120 theta.append((i + 0.5) / (N + 1) * np.pi / 2)
121 theta = np.array(theta)
122 phi = np.mod(phi, 2 * np.pi)
124 if not self._enforce_friedel_symmetry: 124 ↛ 125line 124 didn't jump to line 125, because the condition on line 124 was never true
125 theta = np.concatenate((theta, np.arccos(-np.cos(theta))))
126 phi = np.concatenate((phi, phi))
127 return theta, phi
129 def get_inner_product(self,
130 u: NDArray,
131 v: NDArray) -> NDArray:
132 r""" Retrieves the inner product of two coefficient arrays, that is to say,
133 the sum-product over the last axis.
135 Parameters
136 ----------
137 u
138 The first coefficient array, of arbitrary shape and dimension, so long as
139 the number of coefficients equals the length of this :class:`GaussianKernels` instance.
140 v
141 The second coefficient array, of the same shape as :attr:`u`.
142 """
143 assert u.shape[-1] == len(self)
144 assert u.shape == v.shape
145 return np.einsum('...i, ...i -> ...', u, v, optimize='greedy')
147 def _get_spherical_distances(self,
148 theta_1: NDArray[float], theta_2: NDArray[float],
149 phi_1: NDArray[float], phi_2: NDArray[float],
150 radius: float = 1.,
151 enforce_friedel_symmetry: bool = False) -> NDArray[float]:
152 r""" Function for obtaining the distances between two point sets
153 on a sphere, possibly with Friedel symmetry enforced.
154 Arrays can have any shape, but they must all be broadcastable, and the polar angles of
155 each set must have the same shape as the azimuthal angles.
156 If the first and second set of points have the same shape, then the distances will be
157 computed pointwise. Otherwise, the distances will be computed by broadcasting.
159 Parameters
160 ----------
161 theta_1
162 The polar angle of the first set of points, defined as :math:`\arccos(z_1)`.
163 theta_2
164 The polar angle fo the second set of points, defined as :math:`\arccos(z_2)`.
165 phi_2
166 The azimuthal angle of the first set of points, defined as :math:`\arctan_2(y_1, x_1)`.
167 phi_2
168 The azimuthal angle of the second set of points, defined as :math:`\arctan_2(y_2, x_2)`.
169 radius
170 The radius of the sphere. Default is `1`, i.e., describing a unit sphere.
171 enforce_friedel_symmetry
172 If ``True`` (default), the point :math:`(x, y, z)` will be considered as
173 equivalent to the point :math:`(-x, -y, -z)` and the maximum possible distance on the sphere will
174 be :math:`\frac{\sqrt{x^2 + y^2 + z^2} \pi}{2}`, i.e., at a 90-degree angle.
175 Otherwise, the two points will be considered distinct and the maximum
176 distance will be :math:`\sqrt{x^2 + y^2 + z^2} \pi`.
177 """
178 phi_diff = abs(phi_1 - phi_2)
179 sine_factor = np.sin(theta_1) * np.sin(theta_2) * np.cos(phi_diff)
180 cosine_factor = np.cos(theta_1) * np.cos(theta_2)
181 if enforce_friedel_symmetry: 181 ↛ 184line 181 didn't jump to line 184, because the condition on line 181 was never false
182 return radius * np.arccos(abs(sine_factor + cosine_factor).clip(0., 1.))
183 else:
184 return radius * np.arccos((sine_factor + cosine_factor).clip(-1., 1.))
186 def _get_basis_function_scale_factors(self):
187 """ The basis functions are scaled to have a lower intensity in areas of the
188 half-sphere where the grid is more dense. This function computes those scale factors.
189 """
190 theta, phi = self._get_kurihara_mesh(self._grid_scale)
191 # Probe at location of each kernel function to normalize over sphere
192 mesh_distances = self._get_spherical_distances(
193 theta.reshape(-1, 1), theta.reshape(1, -1),
194 phi.reshape(-1, 1), phi.reshape(1, -1),
195 enforce_friedel_symmetry=self._enforce_friedel_symmetry)
196 std = (self._kernel_scale_parameter * np.sqrt(2 * np.pi)) / (2 * (self._grid_scale + 1))
197 norm_matrix = np.exp(-(1 / 2) * (mesh_distances / std) ** 2)
198 # The normalization factor is the inverse of the unnormalized function value at each grid point
199 norm_factors = np.reciprocal(norm_matrix.sum(-1))
200 return norm_factors
202 def _get_projection_matrix(self, probed_coordinates: ProbedCoordinates = None) -> NDArray[float]:
203 """ Computes the matrix necessary for forward and gradient calculations.
204 Called when the coordinate system has been updated, or one of
205 ``kernel_scale_parameter`` or ``grid_scale`` has been changed."""
206 if probed_coordinates is None: 206 ↛ 207line 206 didn't jump to line 207, because the condition on line 206 was never true
207 probed_coordinates = self.probed_coordinates
208 theta, phi = self._get_kurihara_mesh(self._grid_scale)
209 phi = phi.reshape(1, 1, 1, -1)
210 theta = theta.reshape(1, 1, 1, -1)
211 _, probed_polar_angles, probed_azim_angles = probed_coordinates.to_spherical
212 probed_polar_angles = probed_polar_angles[..., np.newaxis]
213 probed_azim_angles = probed_azim_angles[..., np.newaxis]
214 # Find distances to all probed detector points on sphere
215 distances = self._get_spherical_distances(theta, probed_polar_angles,
216 phi, probed_azim_angles,
217 enforce_friedel_symmetry=self._enforce_friedel_symmetry)
219 # Probe at location of each kernel function to normalize over sphere
220 std = (self._kernel_scale_parameter * np.sqrt(2 * np.pi)) / (2 * (self._grid_scale + 1))
221 matrix = np.exp(-(1 / 2) * (distances / std) ** 2)
223 # The basis functions are scaled to have a lower intensity where the grid is more dense.
224 norm_factors = self._get_basis_function_scale_factors().reshape(1, 1, 1, -1)
225 return matrix * norm_factors
227 def forward(self,
228 coefficients: NDArray,
229 indices: NDArray = None) -> NDArray:
230 """ Carries out a forward computation of projections from Gaussian kernel space
231 into detector space, for one or several tomographic projections.
233 Parameters
234 ----------
235 coefficients
236 An array of coefficients, of arbitrary shape so long as the last
237 axis has the same size as :attr:`~.GaussianKernels.kernel_scale_parameter`, and if
238 :attr:`indices` is ``None`` or greater than one, the first axis should have the
239 same length as :attr:`indices`
240 indices
241 Optional. Indices of the tomographic projections for which the forward
242 computation is to be performed. If ``None``, the forward computation will
243 be performed for all projections.
245 Returns
246 -------
247 An array of values on the detector corresponding to the :attr:`coefficients` given.
248 If :attr:`indices` contains exactly one index, the shape is ``(coefficients.shape[:-1], J)``
249 where ``J`` is the number of detector segments. If :attr:`indices` is ``None`` or contains
250 several indices, the shape is ``(N, coefficients.shape[1:-1], J)`` where ``N``
251 is the number of tomographic projections for which the computation is performed.
253 Notes
254 -----
255 The assumption is made in this implementation that computations over several
256 indices act on sets of images from different projections. For special usage
257 where multiple projections of entire fields are desired, it may be better
258 to use :attr:`projection_matrix` directly. This also applies to
259 :meth:`gradient`.
260 """
261 assert coefficients.shape[-1] == len(self)
262 self._update()
263 output = np.zeros(coefficients.shape[:-1] + (self._projection_matrix.shape[1],),
264 coefficients.dtype)
265 if indices is None:
266 framewise_contraction_transpose(self._projection_matrix,
267 coefficients,
268 output)
269 elif indices.size == 1: 269 ↛ 277line 269 didn't jump to line 277, because the condition on line 269 was never false
270 np.einsum('ijk, ...k -> ...j',
271 self._projection_matrix[indices],
272 coefficients,
273 out=output,
274 optimize='greedy',
275 casting='unsafe')
276 else:
277 framewise_contraction_transpose(self._projection_matrix[indices],
278 coefficients,
279 output)
280 return output
282 def gradient(self,
283 coefficients: NDArray,
284 indices: NDArray = None) -> NDArray:
285 """ Carries out a gradient computation of projections from Gaussian kernel space
286 into detector space for one or several tomographic projections.
288 Parameters
289 ----------
290 coefficients
291 An array of coefficients (or residuals) of arbitrary shape so long as the last
292 axis has the same size as the number of detector segments.
293 indices
294 Optional. Indices of the tomographic projections for which the gradient
295 computation is to be performed. If ``None``, the gradient computation will
296 be performed for all projections.
298 Returns
299 -------
300 An array of gradient values based on the :attr:`coefficients` given.
301 If :attr:`indices` contains exactly one index, the shape is ``(coefficients.shape[:-1], J)``
302 where ``J`` is the number of detector segments. If indices is ``None`` or contains
303 several indices, the shape is ``(N, coefficients.shape[1:-1], J)`` where ``N``
304 is the number of tomographic projections for which the computation is performed.
306 Notes
307 -----
308 When solving an inverse problem, one should not attempt to optimize the
309 coefficients directly using the gradient one obtains by applying this method to the data.
310 Instead, one must either take the gradient of the residual between the
311 :meth:`~.GaussianKernels.forward` computation of the coefficients and the data.
312 Alternatively one can apply both the forward and the gradient computation to the
313 coefficients to be optimized, and the gradient computation to the data, and treat
314 the residual of the two as the gradient of the optimization coefficients. The approaches
315 are algebraically equivalent, but one may be more efficient than the other in some
316 circumstances. However, normally, the projection between detector and
317 ``GaussianKernel`` space is only a small part of the overall computation,
318 so there is typically not much to be gained from optimizing it.
319 """
320 self._update()
321 output = np.zeros(coefficients.shape[:-1] + (self._projection_matrix.shape[2],),
322 coefficients.dtype)
323 if indices is None:
324 framewise_contraction(self._projection_matrix,
325 coefficients,
326 output)
327 elif indices.size == 1: 327 ↛ 335line 327 didn't jump to line 335, because the condition on line 327 was never false
328 np.einsum('ikj, ...k -> ...j',
329 self._projection_matrix[indices],
330 coefficients,
331 out=output,
332 optimize='greedy',
333 casting='unsafe')
334 else:
335 framewise_contraction(self._projection_matrix[indices],
336 coefficients,
337 output)
338 return output
340 def get_amplitudes(self, coefficients: NDArray[float],
341 probed_coordinates: ProbedCoordinates = None) -> NDArray[float]:
342 """ Computes the amplitudes of the spherical function represented by the provided
343 :attr:`coefficients` at the :attr:`probed_coordinates`.
345 Parameters
346 ----------
347 coefficients
348 An array of coefficients of arbitrary shape, provided that the
349 last dimension contains the coefficients for one spherical function.
350 probed_coordinates
351 An instance of :class:`mumott.core.ProbedCoordinates` with its :attr:`vector`
352 attribute indicating the points of the sphere for which to evaluate the amplitudes.
353 """
354 if probed_coordinates is None: 354 ↛ 355line 354 didn't jump to line 355, because the condition on line 354 was never true
355 probed_coordinates = self._probed_coordinates
356 matrix = self._get_projection_matrix(probed_coordinates)
357 self._make_projection_matrix_sparse(matrix)
358 return np.einsum('ij, ...j', matrix.squeeze(), coefficients, optimize='greedy')
360 def get_spherical_harmonic_coefficients(self, coefficients: NDArray, ell_max: int = None):
361 """ Computes the spherical harmonic coefficients of the spherical function
362 represented by the provided :attr:`coefficients` using a Driscoll-Healy grid.
364 For details on the Driscoll-Healy grid, see
365 `the SHTools page <https://shtools.github.io/SHTOOLS/grid-formats.html>`_ for a
366 comprehensive overview.
368 Parameters
369 ----------
370 coefficients
371 An array of coefficients of arbitrary shape, provided that the
372 last dimension contains the coefficients for one function.
373 ell_max
374 The bandlimit of the spherical harmonic expansion. By default, it is ``2 * grid_scale``.
376 """
377 if ell_max is None:
378 ell_max = 2 * self._grid_scale
379 dh_grid_size = 4 * (self._grid_scale + 1)
380 mapper = SphericalHarmonicMapper(ell_max=ell_max, polar_resolution=dh_grid_size,
381 azimuthal_resolution=dh_grid_size,
382 enforce_friedel_symmetry=self._enforce_friedel_symmetry)
383 coordinates = ProbedCoordinates(mapper.unit_vectors.reshape((1, -1, 1, 3)))
384 amplitudes = self.get_amplitudes(coefficients, coordinates)\
385 .reshape(coefficients.shape[:-1] + (dh_grid_size, dh_grid_size))
386 spherical_harmonics_coefficients = mapper.get_harmonic_coefficients(amplitudes)
387 return spherical_harmonics_coefficients
389 def get_second_moments(self, coefficients: NDArray[float]) -> NDArray[float]:
390 """
391 Calculate the second moments of the functions described by :attr:`coefficients`.
393 Parameters
394 ----------
395 coefficients
396 An array of coefficients (or residuals) of arbitrary shape so long as the last
397 axis has the same size as the number of detector channels.
399 Returns
400 -------
401 Array containing the second moments of the functions described by coefficients,
402 formatted as rank-two tensors with tensor indices in the last 2 dimensions.
403 """
405 # Make list of direction vectors
406 theta, phi = self._get_kurihara_mesh(self._grid_scale)
407 direction_vectors = np.stack(
408 (np.sin(theta) * np.cos(phi),
409 np.sin(theta) * np.sin(phi),
410 np.cos(theta),), axis=-1
411 )
413 norm_factors = self._get_basis_function_scale_factors()
414 second_moments_array = np.zeros((*coefficients.shape[:-1], 3, 3))
416 sumint = np.zeros(coefficients.shape[:-1])
417 sumxx = np.zeros(coefficients.shape[:-1])
418 sumxy = np.zeros(coefficients.shape[:-1])
419 sumxz = np.zeros(coefficients.shape[:-1])
420 sumyy = np.zeros(coefficients.shape[:-1])
421 sumyz = np.zeros(coefficients.shape[:-1])
422 sumzz = np.zeros(coefficients.shape[:-1])
424 for mode_number in range(len(self)):
426 sumint += norm_factors[mode_number] * coefficients[..., mode_number]
427 sumxx += norm_factors[mode_number] * coefficients[..., mode_number]\
428 * direction_vectors[mode_number, 0]**2
429 sumxy += norm_factors[mode_number] * coefficients[..., mode_number]\
430 * direction_vectors[mode_number, 0] * direction_vectors[mode_number, 1]
431 sumxz += norm_factors[mode_number] * coefficients[..., mode_number]\
432 * direction_vectors[mode_number, 0] * direction_vectors[mode_number, 2]
433 sumyy += norm_factors[mode_number] * coefficients[..., mode_number]\
434 * direction_vectors[mode_number, 1]**2
435 sumyz += norm_factors[mode_number] * coefficients[..., mode_number]\
436 * direction_vectors[mode_number, 1] * direction_vectors[mode_number, 2]
437 sumzz += norm_factors[mode_number] * coefficients[..., mode_number]\
438 * direction_vectors[mode_number, 2]**2
440 std = (self._kernel_scale_parameter * np.sqrt(2 * np.pi)) / (2 * (self._grid_scale + 1))
441 gaussian_noramlization_term = 1/np.sqrt(2*np.pi*std**2) # This is only approximate on the sphere
443 second_moments_array[..., 0, 0] = sumxx * gaussian_noramlization_term / len(self)
444 second_moments_array[..., 0, 1] = sumxy * gaussian_noramlization_term / len(self)
445 second_moments_array[..., 0, 2] = sumxz * gaussian_noramlization_term / len(self)
446 second_moments_array[..., 1, 0] = sumxy * gaussian_noramlization_term / len(self)
447 second_moments_array[..., 1, 1] = sumyy * gaussian_noramlization_term / len(self)
448 second_moments_array[..., 1, 2] = sumyz * gaussian_noramlization_term / len(self)
449 second_moments_array[..., 2, 0] = sumxz * gaussian_noramlization_term / len(self)
450 second_moments_array[..., 2, 1] = sumyz * gaussian_noramlization_term / len(self)
451 second_moments_array[..., 2, 2] = sumzz * gaussian_noramlization_term / len(self)
453 return second_moments_array
455 def get_output(self,
456 coefficients: NDArray) -> ReconstructionDerivedQuantities:
457 r""" Returns a :class:`ReconstructionDerivedQuantities` instance of output data for
458 a given array of basis set coefficients.
460 Parameters
461 ----------
462 coefficients
463 An array of coefficients of arbitrary shape and dimensions, except
464 its last dimension must be the same length as the :attr:`len` of this instance.
465 Computations only operate over the last axis of :attr:`coefficients`, so derived
466 properties in the output will have the shape ``(*coefficients.shape[:-1], ...)``.
468 Returns
469 -------
470 :class:`ReconstructionDerivedQuantities` containing a number of quantities that
471 have been computed from the spherical functions represented by the input
472 coefficients.
473 """
475 assert coefficients.shape[-1] == len(self)
476 # Update to ensure non-dirty output state.
477 self._update()
479 norm_factors = self._get_basis_function_scale_factors()
480 std = (self._kernel_scale_parameter * np.sqrt(2 * np.pi)) / (2 * (self._grid_scale + 1))
481 gaussian_noramlization_term = 1/np.sqrt(2*np.pi*std**2) # This is only approximate on the sphere
482 mode_integrated_intensities = gaussian_noramlization_term * norm_factors.reshape((1, 1, 1, -1))
484 mean_intensity = np.mean(coefficients * mode_integrated_intensities, axis=-1)
485 second_moment_tensor = self.get_second_moments(coefficients)
486 eigenvalues, eigenvectors = get_sorted_eigenvectors(second_moment_tensor)
487 fractional_anisotropy = np.sqrt((eigenvalues[..., 0] - eigenvalues[..., 1])**2
488 + (eigenvalues[..., 1] - eigenvalues[..., 2])**2
489 + (eigenvalues[..., 2] - eigenvalues[..., 0])**2)
490 fractional_anisotropy = fractional_anisotropy / np.sqrt(2*np.sum(eigenvalues**2, axis=-1))
492 reconstruction_derived_quantities = ReconstructionDerivedQuantities(
493 volume_shape=tuple(coefficients.shape[:3]),
494 mean_intensity=mean_intensity,
495 fractional_anisotropy=fractional_anisotropy,
496 eigenvector_1=np.copy(eigenvectors[..., 0]),
497 eigenvector_2=np.copy(eigenvectors[..., 1]),
498 eigenvector_3=np.copy(eigenvectors[..., 2]),
499 eigenvalue_1=np.copy(eigenvalues[..., 0]),
500 eigenvalue_2=np.copy(eigenvalues[..., 1]),
501 eigenvalue_3=np.copy(eigenvalues[..., 2]),
502 second_moment_tensor=second_moment_tensor
503 )
505 return reconstruction_derived_quantities
507 def __iter__(self) -> Iterator[Tuple[str, Any]]:
508 """ Allows class to be iterated over and in particular be cast as a dictionary.
509 """
510 yield 'name', type(self).__name__
511 yield 'grid_scale', self._grid_scale
512 yield 'kernel_scale_parameter', self._kernel_scale_parameter
513 yield 'enforce_friedel_symmetry', self._enforce_friedel_symmetry
514 yield 'projection_matrix', self._projection_matrix
515 yield 'hash', hex(hash(self))[2:]
517 def __len__(self) -> int:
518 return self._projection_matrix.shape[-1]
520 def __hash__(self) -> int:
521 """Returns a hash reflecting the internal state of the instance.
523 Returns
524 -------
525 A hash of the internal state of the instance,
526 cast as an ``int``.
527 """
528 to_hash = [self._grid_scale,
529 self.grid,
530 self._kernel_scale_parameter,
531 self._enforce_friedel_symmetry,
532 self._projection_matrix,
533 self._probed_coordinates_hash]
534 return int(list_to_hash(to_hash), 16)
536 def _update(self) -> None:
537 # We only run updates if the hashes do not match.
538 if self.is_dirty:
539 self._projection_matrix = self._get_integrated_projection_matrix()
540 self._probed_coordinates_hash = hash(self._probed_coordinates)
542 @property
543 def is_dirty(self) -> bool:
544 return hash(self._probed_coordinates) != self._probed_coordinates_hash
546 @property
547 def projection_matrix(self) -> NDArray:
548 """ The matrix used to project spherical functions from the unit sphere onto the detector.
549 If ``v`` is a vector of gaussian kernel coefficients, and ``M`` is the ``projection_matrix``,
550 then ``M @ v`` gives the corresponding values on the detector segments associated with
551 each projection. ``M[i] @ v`` gives the values on the detector segments associated with
552 projection ``i``.
554 If ``r`` is a residual between a projection from Gaussian kernel to detector space and data from
555 projection ``i``, then ``M[i].T @ r`` gives the associated gradient in Gaussian kernel space.
556 """
557 self._update()
558 return self._projection_matrix
560 @property
561 def grid_scale(self) -> int:
562 """ The number of azimuthal rings from each pole to the equator in the
563 spherical grid.
564 """
565 return self._grid_scale
567 @grid_scale.setter
568 def grid_scale(self, val: int) -> None:
569 if val < 0 or val != round(val):
570 raise ValueError('grid_scale must be a non-negative integer,'
571 f' but a value of {val} was given!')
572 self._grid_scale = val
573 self._projection_matrix = self._get_integrated_projection_matrix()
575 @property
576 def kernel_scale_parameter(self) -> float:
577 """ The scale parameter for each kernel.
578 """
579 return self._kernel_scale_parameter
581 @kernel_scale_parameter.setter
582 def kernel_scale_parameter(self, val: float) -> float:
583 self._kernel_scale_parameter = val
584 self._projection_matrix = self._get_integrated_projection_matrix()
586 @property
587 def enforce_friedel_symmetry(self) -> bool:
588 """ If ``True``, Friedel symmetry is enforced, i.e., the point
589 :math:`-r` is treated as equivalent to :math:`r`. """
590 return self._enforce_friedel_symmetry
592 @property
593 def grid(self) -> Tuple[NDArray['float'], NDArray['float']]:
594 r""" Returns the polar and azimuthal angles of the grid used by the basis.
596 Returns
597 -------
598 A ``Tuple`` with contents ``(polar_angle, azimuthal_angle)``, where the
599 polar angle is defined as :math:`\arccos(z)`.
600 """
601 return self._get_kurihara_mesh(self._grid_scale)
603 @property
604 def grid_hash(self) -> str:
605 """ Returns a hash of :attr:`grid`.
606 """
607 return list_to_hash([self.grid])
609 @property
610 def projection_matrix_hash(self) -> str:
611 """ Returns a hash of :attr:`projection_matrix`.
612 """
613 return list_to_hash([self.projection_matrix])
615 def __str__(self) -> str:
616 wdt = 74
617 s = []
618 s += ['-' * wdt]
619 s += ['GaussianKernels'.center(wdt)]
620 s += ['-' * wdt]
621 with np.printoptions(threshold=4, edgeitems=2, precision=5, linewidth=60):
622 s += ['{:18} : {}'.format('grid_scale', self.grid_scale)]
623 s += ['{:18} : {}'.format('grid_hash', self.grid_hash)]
624 s += ['{:18} : {}'.format('enforce_friedel_symmetry', self.enforce_friedel_symmetry)]
625 s += ['{:18} : {}'.format('kernel_scale_parameter', self.kernel_scale_parameter)]
626 s += ['{:18} : {}'.format('projection_matrix_hash', self.projection_matrix_hash)]
627 s += ['{:18} : {}'.format('hash', hex(hash(self))[2:8])]
628 s += ['-' * wdt]
629 return '\n'.join(s)
631 def _repr_html_(self) -> str:
632 s = []
633 s += [f'<h3>{self.__class__.__name__}</h3>']
634 s += ['<table border="1" class="dataframe">']
635 s += ['<thead><tr><th style="text-align: left;">Field</th><th>Size</th><th>Data</th></tr></thead>']
636 s += ['<tbody>']
637 with np.printoptions(threshold=4, edgeitems=2, precision=2, linewidth=40):
638 s += ['<tr><td style="text-align: left;">grid_scale</td>']
639 s += [f'<td>1</td><td>{self.grid_scale}</td></tr>']
640 s += ['<tr><td style="text-align: left;">grid_hash</td>']
641 s += [f'<td>{len(self.grid_hash)}</td><td>{self.grid_hash[:6]}</td></tr>']
642 s += ['<tr><td style="text-align: left;">kernel_scale_parameter</td>']
643 s += [f'<td>1</td><td>{self.kernel_scale_parameter}</td></tr>']
644 s += ['<tr><td style="text-align: left;">enforce_friedel_symmetry</td>']
645 s += [f'<td>1</td>'
646 f'<td>{self.enforce_friedel_symmetry}</td></tr>']
647 s += ['<tr><td style="text-align: left;">projection_matrix</td>']
648 s += [f'<td>{len(self.projection_matrix_hash)}</td>'
649 f'<td>{self.projection_matrix_hash[:6]}</td></tr>']
650 s += ['<tr><td style="text-align: left;">hash</td>']
651 s += [f'<td>{len(hex(hash(self)))}</td><td>{hex(hash(self))[2:8]}</td></tr>']
652 s += ['</tbody>']
653 s += ['</table>']
654 return '\n'.join(s)