Coverage for local_installation_linux/mumott/methods/basis_sets/spherical_harmonics.py: 93%
190 statements
« prev ^ index » next coverage.py v7.3.2, created at 2025-05-05 21:21 +0000
1import logging
2from typing import Any, Iterator, List, Tuple
3from numpy.typing import NDArray
5import numpy as np
6from scipy.special import sph_harm
8from mumott import ProbedCoordinates
9from mumott.core.hashing import list_to_hash
10from mumott.methods.utilities.tensor_operations import (framewise_contraction,
11 framewise_contraction_transpose)
12from mumott.output_handling.reconstruction_derived_quantities import\
13 (ReconstructionDerivedQuantities, get_sorted_eigenvectors)
14from .base_basis_set import BasisSet
17logger = logging.getLogger(__name__)
20class SphericalHarmonics(BasisSet):
21 """ Basis set class for spherical harmonics, the canonical representation
22 of polynomials on the unit sphere and a simple way of representing
23 band-limited spherical functions which allows for easy computations of statistics
24 and is suitable for analyzing certain symmetries.
26 Parameters
27 ----------
28 ell_max : int
29 The bandlimit of the spherical functions that you want to be able to represent.
30 A good rule of thumb is that :attr:`ell_max` should not exceed the
31 number of detector segments minus 1.
32 probed_coordinates : ProbedCoordinates
33 Optional. A container with the coordinates on the sphere probed at each detector segment by the
34 experimental method. Its construction from the system geometry is method-dependent.
35 By default, an empty instance of
36 :class:`ProbedCoordinates <mumott.core.probed_coordinates.ProbedCoordinates>` is created.
37 enforce_friedel_symmetry : bool
38 If set to ``True``, Friedel symmetry will be enforced, using the assumption that points
39 on opposite sides of the sphere are equivalent. This results in only even ``ell`` being used.
40 kwargs
41 Miscellaneous arguments which relate to segment integrations can be
42 passed as keyword arguments:
44 integration_mode
45 Mode to integrate line segments on the reciprocal space sphere. Possible options are
46 ``'simpson'``, ``'midpoint'``, ``'romberg'``, ``'trapezoid'``.
47 ``'simpson'``, ``'trapezoid'``, and ``'romberg'`` use adaptive
48 integration with the respective quadrature rule from ``scipy.integrate``.
49 ``'midpoint'`` uses a single mid-point approximation of the integral.
50 Default value is ``'simpson'``.
51 n_integration_starting_points
52 Number of points used in the first iteration of the adaptive integration.
53 The number increases by the rule ``N`` ← ``2 * N - 1`` for each iteration.
54 Default value is 3.
55 integration_tolerance
56 Tolerance for the maximum relative error between iterations before the integral
57 is considered converged. Default is ``1e-5``.
58 integration_maxiter
59 Maximum number of iterations. Default is ``10``.
60 """
61 def __init__(self,
62 probed_coordinates: ProbedCoordinates = None,
63 ell_max: int = 0,
64 enforce_friedel_symmetry: bool = True,
65 **kwargs):
66 super().__init__(probed_coordinates, **kwargs)
67 self._probed_coordinates_hash = hash(self.probed_coordinates)
68 self._ell_max = ell_max
69 self._ell_indices = np.zeros(1)
70 self._emm_indices = np.zeros(1)
71 # Compute initial values for indices and matrix.
72 self._enforce_friedel_symmetry = enforce_friedel_symmetry
73 self._calculate_coefficient_indices()
74 self._projection_matrix = self._get_integrated_projection_matrix()
76 def _calculate_coefficient_indices(self) -> None:
77 """
78 Computes the attributes :attr:`~.SphericalHarmonics.ell_indices` and
79 :attr:`~.SphericalHarmonics.emm_indices`. Called when :attr:`~.SphericalHarmonics.ell_max`
80 changes.
81 """
82 if self._enforce_friedel_symmetry:
83 divisor = 2
84 else:
85 divisor = 1
87 mm = np.zeros((self._ell_max + 1) * (self._ell_max // divisor + 1), dtype=int)
88 ll = np.zeros((self._ell_max + 1) * (self._ell_max // divisor + 1), dtype=int)
89 count = 0
90 for h in range(0, self._ell_max + 1, divisor):
91 for i in range(-h, h + 1):
92 ll[count] = h
93 mm[count] = i
94 count += 1
95 self._ell_indices = ll
96 self._emm_indices = mm
98 def _get_projection_matrix(self, probed_coordinates: ProbedCoordinates = None) -> None:
99 """ Computes the matrix necessary for forward and gradient calculations.
100 Called when the coordinate system has been updated or ``ell_max`` has changed."""
101 if probed_coordinates is None: 101 ↛ 102line 101 didn't jump to line 102, because the condition on line 101 was never true
102 probed_coordinates = self._probed_coordinates
103 _, probed_polar_angles, probed_azim_angles = probed_coordinates.to_spherical
104 # retrieve complex spherical harmonics with emm >= 0, shape (N, M, len(ell_indices), K)
105 complex_factors = sph_harm(abs(self._emm_indices)[np.newaxis, np.newaxis, np.newaxis, ...],
106 self._ell_indices[np.newaxis, np.newaxis, np.newaxis, ...],
107 probed_azim_angles[..., np.newaxis],
108 probed_polar_angles[..., np.newaxis])
109 # cancel Condon-Shortley phase factor in scipy.special.sph_harm
110 condon_shortley_factor = (-1.) ** self._emm_indices
111 # 4pi normalization factor and complex-to-real normalization factor for m != 0
112 norm_factor = np.sqrt(4 * np.pi) * \
113 np.sqrt(1 + (self._emm_indices != 0).astype(int)) * condon_shortley_factor
114 matrix = norm_factor[np.newaxis, np.newaxis, np.newaxis, ...] * (
115 (self._emm_indices >= 0)[np.newaxis, np.newaxis, np.newaxis, ...] * complex_factors.real +
116 (self._emm_indices < 0)[np.newaxis, np.newaxis, np.newaxis, ...] * complex_factors.imag)
117 return matrix
119 def forward(self,
120 coefficients: np.array,
121 indices: np.array = None) -> np.array:
122 """ Carries out a forward computation of projections from spherical harmonic space
123 into detector space, for one or several tomographic projections.
125 Parameters
126 ----------
127 coefficients
128 An array of coefficients, of arbitrary shape so long as the last
129 axis has the same size as :attr:`~.SphericalHarmonics.ell_indices`, and if
130 :attr:`indices` is `None` or greater than one, the first axis should have the
131 same length as :attr:`indices`
132 indices
133 Optional. Indices of the tomographic projections for which the forward
134 computation is to be performed. If ``None``, the forward computation will
135 be performed for all projections.
137 Returns
138 -------
139 An array of values on the detector corresponding to the :attr:`coefficients` given.
140 If :attr:`indices` contains exactly one index, the shape is ``(coefficients.shape[:-1], J)``
141 where ``J`` is the number of detector segments. If :attr:`indices` is ``None`` or contains
142 several indices, the shape is ``(N, coefficients.shape[1:-1], J)`` where ``N``
143 is the number of tomographic projections for which the computation is performed.
145 Notes
146 -----
147 The assumption is made in this implementation that computations over several
148 indices act on sets of images from different projections. For special usage
149 where multiple projections of entire fields is desired, it may be better
150 to use :attr:`projection_matrix` directly. This also applies to
151 :meth:`gradient`.
152 """
153 assert coefficients.shape[-1] == self._ell_indices.size
154 self._update()
155 output = np.zeros(coefficients.shape[:-1] + (self._projection_matrix.shape[1],),
156 coefficients.dtype)
157 if indices is None:
158 framewise_contraction_transpose(self._projection_matrix,
159 coefficients,
160 output)
161 elif indices.size == 1: 161 ↛ 169line 161 didn't jump to line 169, because the condition on line 161 was never false
162 np.einsum('ijk, ...k -> ...j',
163 self._projection_matrix[indices],
164 coefficients,
165 out=output,
166 optimize='greedy',
167 casting='unsafe')
168 else:
169 framewise_contraction_transpose(self._projection_matrix[indices],
170 coefficients,
171 output)
172 return output
174 def gradient(self,
175 coefficients: np.array,
176 indices: np.array = None) -> np.array:
177 """ Carries out a gradient computation of projections from spherical harmonic space
178 into detector space, for one or several tomographic projections.
180 Parameters
181 ----------
182 coefficients
183 An array of coefficients (or residuals) of arbitrary shape so long as the last
184 axis has the same size as the number of detector segments.
185 indices
186 Optional. Indices of the tomographic projections for which the gradient
187 computation is to be performed. If ``None``, the gradient computation will
188 be performed for all projections.
190 Returns
191 -------
192 An array of gradient values based on the `coefficients` given.
193 If :attr:`indices` contains exactly one index, the shape is ``(coefficients.shape[:-1], J)``
194 where ``J`` is the number of detector segments. If indices is ``None`` or contains
195 several indices, the shape is ``(N, coefficients.shape[1:-1], J)`` where ``N``
196 is the number of tomographic projections for which the computation is performed.
198 Notes
199 -----
200 When solving an inverse problem, one should not to attempt to optimize the
201 coefficients directly using the ``gradient`` one obtains by applying this method to the data.
202 Instead, one must either take the gradient of the residual between the
203 :meth:`~.SphericalHarmonics.forward` computation of the coefficients and the data.
204 Alternatively one can apply both the forward and the gradient computation to the
205 coefficients to be optimized, and the gradient computation to the data, and treat
206 the residual of the two as the gradient of the optimization coefficients. The approaches
207 are algebraically equivalent, but one may be more efficient than the other in some
208 circumstances.
209 """
210 self._update()
211 output = np.zeros(coefficients.shape[:-1] + (self._projection_matrix.shape[2],),
212 coefficients.dtype)
213 if indices is None:
214 framewise_contraction(self._projection_matrix,
215 coefficients,
216 output)
217 elif indices.size == 1: 217 ↛ 225line 217 didn't jump to line 225, because the condition on line 217 was never false
218 np.einsum('ikj, ...k -> ...j',
219 self._projection_matrix[indices],
220 coefficients,
221 out=output,
222 optimize='greedy',
223 casting='unsafe')
224 else:
225 framewise_contraction(self._projection_matrix[indices],
226 coefficients,
227 output)
228 return output
230 def get_inner_product(self,
231 u: np.array,
232 v: np.array,
233 resolve_spectrum: bool = False,
234 spectral_moments: List[int] = None) -> np.array:
235 r""" Retrieves the inner product of two coefficient arrays.
237 Notes
238 -----
239 The canonical inner product in a spherical harmonic representation
240 is :math:`\sum_\ell N(\ell) \sum_m u_m^\ell v_m^\ell`, where :math:`N(\ell)` is
241 a normalization constant (which is unity for the :math:`4\pi` normalization).
242 This inner product is a rotational invariant. The rotational invariance also holds
243 for any partial sums over :math:`\ell`.
244 One can define a function of :math:`\ell` that returns such
245 products, namely :math:`S(\ell, u, v) = N(\ell)\sum_m u_m^\ell v_m^\ell`,
246 called the spectral power function.
247 The sum :math:`\sum_{\ell = 1}S(\ell)` is equal to the covariance of the
248 band-limited spherical functions represented by :math:`u` and :math:`v`, and each
249 :math:`S(\ell, u, v)` is the contribution to the covariance of the band :math:`\ell`.
250 See also
251 `the SHTOOLS documentation <https://shtools.github.io/SHTOOLS/real-spherical-harmonics.html>`_
252 for an excellent overview of this.
254 Parameters
255 ----------
256 u
257 The first coefficient array, of arbitrary shape and dimension,
258 except the last dimension must be the same as the length of
259 :attr:`~.SphericalHarmonics.ell_indices`.
260 v
261 The second coefficient array, of the same shape as ``u``.
262 resolve_spectrum
263 Optional. Whether to resolve the product according to each frequency band, given by the
264 coefficients of each ``ell`` in :attr:`~SphericalHarmonics.ell_indices`.
265 Defaults to ``False``, which means that the sum of every component of the spectrum is returned.
266 If ``True``, components are returned in order of ascending ``ell``.
267 The ``ell`` included in the spectrum depends on :attr:`spectral_moments`.
268 spectral_moments
269 Optional. List of particular values of ``ell`` to calculate the inner product for.
270 Defaults to ``None``, which is identical to including all values of ``ell``
271 in the calculation.
272 If :attr:`spectral_moments` contains all nonzero values of ``ell`` and
273 :attr:`resolve_spectrum` is ``False``, the covariance of
274 :attr:`v` and :attr:`u` will be calculated (the sum of the inner product over all non-zero ``ell``
275 If ``resolve_spectrum`` is ``True``, the covariance per `ell` in ``spectral_moments``, will
276 be calculated, i.e., the inner products will not be summed over.
278 Returns
279 -------
280 An array of the inner products of the spherical functions represented by ``u`` and ``v``.
281 Has the shape ``(u.shape[:-1])`` if :attr:`resolve_spectrum` is ``False``,
282 ``(u.shape[:-1] + (ell_max // 2 + 1,))`` if :attr:`resolve_spectrum` is ``True`` and
283 ``spectral_moments`` is ``None``, and finally the shape
284 ``(u.shape[:-1] + (np.unique(spectral_moments).size,))`` if :attr:`resolve_spectrum` is ``True``
285 and :attr:`spectral_moments` is a list of integers found
286 in :attr:`~.SphericalHarmonics.ell_indices`
287 """
288 assert u.shape == v.shape
289 assert u.shape[-1] == self._ell_indices.size
290 if not resolve_spectrum:
291 if spectral_moments is None:
292 return np.einsum('...i, ...i -> ...', u, v)
293 # pick out only the subset where ell matches the provided spectral moments
294 where = np.any([np.equal(self._ell_indices, ell) for ell in spectral_moments], axis=0)
295 return np.einsum('...i, ...i -> ...', u[..., where], v[..., where])
296 if spectral_moments is None: 296 ↛ 299line 296 didn't jump to line 299, because the condition on line 296 was never false
297 which_ell = np.unique(self._ell_indices)
298 else:
299 which_ell = np.unique(spectral_moments)
300 which_ell = [ell for ell in which_ell if np.any(np.equals(self._ell_indices, ell))]
301 power_spectrum = np.zeros((*u.shape[:-1], which_ell.size))
302 # power spectrum for any one ell is given by inner product over each ell
303 for i, ell in enumerate(which_ell):
304 power_spectrum[..., i] = np.einsum('...i, ...i -> ...',
305 u[..., self._ell_indices == ell],
306 v[..., self._ell_indices == ell],
307 optimize='greedy')
308 return power_spectrum
310 def get_covariances(self,
311 u: np.array,
312 v: np.array,
313 resolve_spectrum: bool = False) -> np.array:
314 """ Returns the covariances of the spherical functions represented by
315 two coefficient arrays.
317 Parameters
318 ----------
319 u
320 The first coefficient array, of arbitrary shape
321 except its last dimension must be the same length as
322 the length of :attr:`~SphericalHarmonics.ell_indices`.
323 v
324 The second coefficient array, of the same shape as :attr:`u`.
325 resolve_spectrum
326 Optional. Whether to resolve the product according to each frequency band, given by the
327 coefficients of each ``ell`` in :attr:`~.SphericalHarmonics.ell_indices`.
328 Default value is ``False``.
330 Returns
331 -------
332 An array of the covariances of the spherical functions represented by ``u`` and ``v``.
333 Has the shape ``(u.shape[:-1])`` if `resolve_spectrum` is ``False``, and
334 ``(u.shape[:-1] + (ell_max // 2 + 1,))`` if `resolve_spectrum` is ``True``, where
335 ``ell_max`` is :attr:`.SphericalHarmonics.ell_max`.
337 Notes
338 -----
339 Calling this function is equivalent to calling :func:`~.SphericalHarmonics.get_inner_product`
340 with ``spectral_moments=np.unique(ell_indices[ell_indices > 0])`` where ``ell_indices`` is
341 :attr:`.SphericalHarmonics.ell_indices`. See the note to
342 :func:`~.SphericalHarmonics.get_inner_product` for mathematical details.
343 """
344 spectral_moments = np.unique(self._ell_indices[self._ell_indices > 0])
345 return self.get_inner_products(u, v, resolve_spectrum, spectral_moments)
347 def get_output(self,
348 coefficients: np.array) -> ReconstructionDerivedQuantities:
349 r""" Returns a :class:`ReconstructionDerivedQuantities` instance of output data for
350 a given array of basis set coefficients.
352 Parameters
353 ----------
354 coefficients
355 An array of coefficients of arbitrary shape and dimensions, except
356 its last dimension must be the same length as the :attr:`len` of this instance.
357 Computations only operate over the last axis of :attr:`coefficients`, so derived
358 properties in the output will have the shape ``(*coefficients.shape[:-1], ...)``.
360 Returns
361 -------
362 :class:`ReconstructionDerivedQuantities` containing a number of quantities that
363 have been computed from the spherical functions represented by the input
364 coefficients.
365 """
367 assert coefficients.shape[-1] == len(self)
368 # Update to ensure non-dirty output state.
369 self._update()
371 mean_intensity = coefficients[..., 0]
372 second_moment_tensor = self._get_rank2_tensor(coefficients)
373 eigenvalues, eigenvectors = get_sorted_eigenvectors(second_moment_tensor)
374 fractional_anisotropy = np.sqrt((eigenvalues[..., 0] - eigenvalues[..., 1])**2
375 + (eigenvalues[..., 1] - eigenvalues[..., 2])**2
376 + (eigenvalues[..., 2] - eigenvalues[..., 0])**2)
377 fractional_anisotropy = fractional_anisotropy / np.sqrt(2*np.sum(eigenvalues**2, axis=-1))
379 reconstruction_derived_quantities = ReconstructionDerivedQuantities(
380 volume_shape=tuple(coefficients.shape[:3]),
381 mean_intensity=mean_intensity,
382 fractional_anisotropy=fractional_anisotropy,
383 eigenvector_1=np.copy(eigenvectors[..., 0]),
384 eigenvector_2=np.copy(eigenvectors[..., 1]),
385 eigenvector_3=np.copy(eigenvectors[..., 2]),
386 eigenvalue_1=np.copy(eigenvalues[..., 0]),
387 eigenvalue_2=np.copy(eigenvalues[..., 1]),
388 eigenvalue_3=np.copy(eigenvalues[..., 2]),
389 second_moment_tensor=second_moment_tensor
390 )
392 return reconstruction_derived_quantities
394 def get_spherical_harmonic_coefficients(
395 self,
396 coefficients: np.array,
397 ell_max: int = None) -> np.array:
398 """ Convert a set of spherical harmonics coefficients to a different :attr:`ell_max`
399 by either zero-padding or truncation and return the result.
401 Parameters
402 ----------
403 coefficients
404 An array of coefficients of arbitrary shape, provided that the
405 last dimension contains the coefficients for one function.
406 ell_max
407 The band limit of the spherical harmonic expansion.
408 """
410 if coefficients.shape[-1] != len(self): 410 ↛ 411line 410 didn't jump to line 411, because the condition on line 410 was never true
411 raise ValueError(f'The number of coefficients ({coefficients.shape[-1]}) does not '
412 f'match the expected value. ({len(self)})')
414 if self._enforce_friedel_symmetry:
415 num_coeff_output = (ell_max+1) * (ell_max+2) // 2
416 elif not self._enforce_friedel_symmetry: 416 ↛ 419line 416 didn't jump to line 419, because the condition on line 416 was never false
417 num_coeff_output = (ell_max+1)**2
419 output_array = np.zeros((*coefficients.shape[:-1], num_coeff_output))
420 output_array[..., :min(len(self), num_coeff_output)] = \
421 coefficients[..., :min(len(self), num_coeff_output)]
422 return output_array
424 def _get_rank2_tensor(self, coefficients: np.array(float)) -> np.array(float):
425 """ Computes the second moments of the reciprocal space maps
426 represented by `coefficients`.
428 Parameters
429 ----------
430 coefficients
431 An array of coefficients of arbitrary shape, as long as the last axis has the same shape
432 as :attr:`~.SphericalHarmonics.ell_indices`.
434 Returns
435 -------
436 First, the 2nd moment tensors of all the functions represented by the :attr:`coefficients`,
437 in 3-by-3 matrix form stored in the last two indices:
438 ``[[xx, xy, xz], [xy, yy, yz], [xz, yz, zz]]``
439 """
440 r2_tensor = coefficients[..., 0, np.newaxis, np.newaxis]\
441 * np.eye(3)[np.newaxis, np.newaxis, np.newaxis, :, :] * 1/3
443 if self._ell_max < 2:
444 logger.info('Note: ell_max < 2, so rank-2 tensors will all be diagonal matrices.')
445 return r2_tensor
447 r2_tensor[..., 0, 1] = coefficients[..., 1] / np.sqrt(15)
448 r2_tensor[..., 1, 0] = coefficients[..., 1] / np.sqrt(15)
449 r2_tensor[..., 1, 2] = coefficients[..., 2] / np.sqrt(15)
450 r2_tensor[..., 2, 1] = coefficients[..., 2] / np.sqrt(15)
451 r2_tensor[..., 2, 2] += coefficients[..., 3] * 2 * np.sqrt(5) / 15
452 r2_tensor[..., 0, 0] += -coefficients[..., 3] * np.sqrt(5) / 15 + coefficients[..., 5] / np.sqrt(15)
453 r2_tensor[..., 1, 1] += -coefficients[..., 3] * np.sqrt(5) / 15 - coefficients[..., 5] / np.sqrt(15)
454 r2_tensor[..., 0, 2] = coefficients[..., 4] / np.sqrt(15)
455 r2_tensor[..., 2, 0] = coefficients[..., 4] / np.sqrt(15)
457 return r2_tensor
459 def __iter__(self) -> Iterator[Tuple[str, Any]]:
460 """ Allows class to be iterated over and in particular be cast as a dictionary.
461 """
462 yield 'name', type(self).__name__
463 yield 'ell_max', self._ell_max
464 yield 'ell_indices', self._ell_indices
465 yield 'emm_indices', self._emm_indices
466 yield 'projection_matrix', self._projection_matrix
467 yield 'hash', hex(hash(self))[2:]
469 def __len__(self) -> int:
470 return len(self._ell_indices)
472 def __hash__(self) -> int:
473 """Returns a hash reflecting the internal state of the instance.
475 Returns
476 -------
477 A hash of the internal state of the instance,
478 cast as an ``int``.
479 """
480 to_hash = [self._ell_max,
481 self._ell_indices,
482 self._emm_indices,
483 self._projection_matrix,
484 self._probed_coordinates_hash]
485 return int(list_to_hash(to_hash), 16)
487 def _update(self) -> None:
488 # We only run updates if the hashes do not match.
489 if self.is_dirty:
490 self._projection_matrix = self._get_integrated_projection_matrix()
491 self._probed_coordinates_hash = hash(self._probed_coordinates)
493 @property
494 def is_dirty(self) -> bool:
495 return hash(self._probed_coordinates) != self._probed_coordinates_hash
497 @property
498 def projection_matrix(self) -> np.array:
499 """ The matrix used to project spherical functions from the unit sphere onto the detector.
500 If ``v`` is a vector of spherical harmonic coefficients, and ``M`` is the ``projection_matrix``,
501 then ``M @ v`` gives the corresponding values on the detector segments associated with
502 each projection. ``M[i] @ v`` gives the values on the detector segments associated with
503 projection ``i``.
505 If ``r`` is a residual between a projection from spherical to detector space and data from
506 projection ``i``, then ``M[i].T @ r`` gives the associated gradient in spherical harmonic
507 space.
508 """
509 self._update()
510 return self._projection_matrix
512 @property
513 def ell_max(self) -> int:
514 r""" The maximum ``ell`` used to represent spherical functions.
516 Notes
517 -----
518 The word ``ell`` is used to represent the cursive small L, also written
519 :math:`\ell`, often used as an index for the degree of the Legendre polynomial
520 in the definition of the spherical harmonics.
521 """
522 return self._ell_max
524 @ell_max.setter
525 def ell_max(self, val: int) -> np.array(int):
526 if (val % 2 != 0 and self._enforce_friedel_symmetry) or val < 0 or val != round(val):
527 raise ValueError('ell_max must be an even (if Friedel symmetry is enforced),'
528 ' non-negative integer,'
529 f' but a value of {val} was given!')
530 self._ell_max = val
531 self._calculate_coefficient_indices()
532 self._projection_matrix = self._get_integrated_projection_matrix()
534 @property
535 def ell_indices(self) -> NDArray[int]:
536 r""" The ``ell`` associated with each coefficient and its corresponding
537 spherical harmonic. Updated when :attr:`~.SphericalHarmonics.ell_max` changes.
539 Notes
540 -----
541 The word ``ell`` is used to represent the cursive small L, also written
542 :math:`\ell`, often used as an index for the degree of the Legendre polynomial
543 in the definition of the spherical harmonics.
544 """
545 return self._ell_indices
547 @property
548 def emm_indices(self) -> NDArray[int]:
549 r""" The ``emm`` associated with each coefficient and its corresponding
550 spherical harmonic. Updated when :attr:`~.SphericalHarmonics.ell_max` changes.
552 Notes
553 -----
554 For consistency with :attr:`~.SphericalHarmonics.ell_indices`, and to avoid
555 visual confusion with other letters, ``emm`` is
556 used to represent the index commonly written :math:`m` in mathematical notation,
557 the frequency of the sine-cosine parts of the spherical harmonics,
558 often called the spherical harmonic order.
559 """
560 return self._emm_indices
562 def __str__(self) -> str:
563 wdt = 74
564 s = []
565 s += ['-' * wdt]
566 s += ['SphericalHarmonics'.center(wdt)]
567 s += ['-' * wdt]
568 with np.printoptions(threshold=4, edgeitems=2, precision=5, linewidth=60):
569 s += ['{:18} : {}'.format('Maximum "ell"', self.ell_max)]
570 s += ['{:18} : {}'.format('"ell" indices', self.ell_indices)]
571 s += ['{:18} : {}'.format('"emm" indices', self.emm_indices)]
572 s += ['{:18} : {}'.format('Projection matrix', self.projection_matrix)]
573 s += ['{:18} : {}'.format('Hash', hex(hash(self))[2:8])]
574 s += ['-' * wdt]
575 return '\n'.join(s)
577 def _repr_html_(self) -> str:
578 s = []
579 s += ['<h3>SphericalHarmonics</h3>']
580 s += ['<table border="1" class="dataframe">']
581 s += ['<thead><tr><th style="text-align: left;">Field</th><th>Size</th><th>Data</th></tr></thead>']
582 s += ['<tbody>']
583 with np.printoptions(threshold=4, edgeitems=2, precision=2, linewidth=40):
584 s += ['<tr><td style="text-align: left;">Maximum "ell"</td>']
585 s += [f'<td>{1}</td><td>{self.ell_max}</td></tr>']
586 s += ['<tr><td style="text-align: left;">"ell" indices</td>']
587 s += [f'<td>{len(self.ell_indices)}</td><td>{self.ell_indices}</td></tr>']
588 s += ['<tr><td style="text-align: left;">"emm" indices</td>']
589 s += [f'<td>{len(self.emm_indices)}</td><td>{self.emm_indices}</td></tr>']
590 s += ['<tr><td style="text-align: left;">Coefficient projection matrix</td>']
591 s += [f'<td>{self.projection_matrix.shape}</td><td>{self.projection_matrix}</td></tr>']
592 s += ['<tr><td style="text-align: left;">Hash</td>']
593 s += [f'<td>{len(hex(hash(self)))}</td><td>{hex(hash(self))[2:8]}</td></tr>']
594 s += ['</tbody>']
595 s += ['</table>']
596 return '\n'.join(s)