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