Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement different normalization and phase conventions for spherical harmonic basis matrices #102

Open
wants to merge 15 commits into
base: develop
Choose a base branch
from

Conversation

mberz
Copy link
Member

@mberz mberz commented Feb 20, 2025

Taken from #82

@mberz mberz added the enhancement New feature or request label Feb 21, 2025
@mberz mberz added this to the v1.0.0 milestone Feb 21, 2025
@tluebeck tluebeck requested review from f-brinkmann and ahms5 March 5, 2025 20:48
Copy link
Member Author

@mberz mberz left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for cleaning this up and implementing the tests.
I can only comment since I opened the PR ;)

I have a couple of comments. Some I only noted in single places, but apply to all functions. It's probably easy to spot.

Could you especially please add a test against the old reference files to ensure that the new "default" settings still agree with the old behavior of the functions.

Comment on lines 389 to 395
if not isinstance(coordinates, pf.Coordinates):
axis = np.where(coordinates.shape == 3)[0][0]
if axis == 0:
coordinates = coordinates.T
coordinates = pf.Coordinates(coordinates[:, 0],
coordinates[:, 1],
coordinates[:, 2])
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think we should do this transform. We should stick to allowing only pf.Coordinates or Samplings here.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It was a leftover from the old branch which I did not touch. But I removed it now.

Comment on lines +317 to +332
if channel_convention == "fuma":
order, degree = fuma_to_nm(acn)
else:
order, degree = acn_to_nm(acn)
basis[:, acn] = _special.spherical_harmonic(
order,
degree,
coords.colatitude,
coords.azimuth)

order, degree, coordinates.colatitude, coordinates.azimuth
)
if normalization == "sn3d":
basis[:, acn] *= n3d_to_sn3d_norm(order)
elif normalization == "maxN":
basis[:, acn] *= n3d_to_maxn(acn)
if phase_convention is None:
# Condon-Shortley phase term is already included in
# the special.spherical_harmonic function
# so need to divide by (-1)^m
basis[:, acn] /= (-1) ** float(degree)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Couldn't these re-normalization and re-sorting operations be moved out of the loop? This should be much faster.

Copy link
Member

@f-brinkmann f-brinkmann Mar 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

+1 for this. I think scipy.special.sph_har can also generate the entire matrix without being called in the loop, but that would require refactoring our special module. Could be a speed-up for high orders and large coordinates objects

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I can see the benefit but the gains are minimal if you are looking at spherical harmonics of order < 10 as is typical in acoustics. I don't think its necessary to change imho.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With the same logic we could make this "jit"-able, but it's a bit of an overkill.

Comment on lines 464 to 479
if channel_convention == "fuma":
order, degree = fuma_to_nm(acn)
else:
order, degree = acn_to_nm(acn)
basis[:, acn] = _special.spherical_harmonic_real(
order,
degree,
coords.colatitude,
coords.azimuth)
order, degree, coordinates.colatitude, coordinates.azimuth
)
if normalization == "sn3d":
basis[:, acn] *= n3d_to_sn3d_norm(order)
elif normalization == "maxN":
basis[:, acn] *= n3d_to_maxn(acn)
if phase_convention is None:
# Condon-Shortley phase term is already included in
# the special.spherical_harmonic function
# so need to divide by (-1)^m
basis[:, acn] /= (-1) ** float(degree)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same comment as above, if possible, please move the re-normalization and sorting out of the loop.

Comment on lines 475 to 479
if phase_convention is None:
# Condon-Shortley phase term is already included in
# the special.spherical_harmonic function
# so need to divide by (-1)^m
basis[:, acn] /= (-1) ** float(degree)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the spherical_harmonic_real function does not include the CS phase. Please check the issue I created on this and compare.
Technically, there's no difference between multiplying with (-1)^m or dividing by, but the comment is not correct then ;)

Comment on lines 534 to 540
if not isinstance(coordinates, pf.Coordinates):
axis = np.where(coordinates.shape == 3)[0][0]
if axis == 0:
coordinates = coordinates.T
coordinates = pf.Coordinates(coordinates[:, 0],
coordinates[:, 1],
coordinates[:, 2])
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same comment as above. I don't think we should do auto-transforming coordinates.,

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as above. It is removed now.

Comment on lines 73 to 81
Y = np.genfromtxt(f'./tests/data/Y_real_{phase_convention}_'
f'{normalization}_{channel_convention}.csv',
dtype=float,
delimiter=',')
basis = sh.spherical_harmonic_basis_real(n_max, coords,
normalization,
channel_convention,
phase_convention=phase_convention)
np.testing.assert_allclose(basis, Y, atol=1e-13)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you also for the benefit of my sanity keep a test against the old files. Just so we're sure that we don't accidentally change the default normaltion and phase conventions of the files? ;)

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sure, I added test against your old n10 files for both, real and complex. However, the new default for Condon-Shortley is None, so we explicitly need to test against phase_convention='Condon-Shortley', not the default.. Probably we should discuss at this point if we want to have "Condon-Shortley" as the default, to be consistent with old implementations?

return basis


def spherical_harmonic_basis_gradient(n_max, coords):
def spherical_harmonic_basis_gradient(n_max, coordinates):
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You did not implement the phase and normalization conventions for the gradients. Is there a reason for this?
Did you want to keep this for a separate PR?

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No, just not thought about it :D. I will implement it in this branch.

normalization : str, optional
Normalization convention, either 'n3d', 'maxN' or 'sn3d'.
The default is 'n3d'.
(maxN is only supported up to 3rd order)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since this is only valid up to third order we should raise an error on invalid inputs.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

channel_convention : str, optional
Channel ordering convention, either 'acn' or 'fuma'.
The default is 'acn'.
(FuMa is only supported up to 3rd order)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since this is only valid up to third order we should raise an error on invalid inputs.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

@mberz mberz requested a review from xefonon March 14, 2025 13:34

Y_n^m(\theta, \phi) = \sqrt{\frac{2n+1}{4\pi}
\frac{(n-m)!}{(n+m)!}} P_n^m(\cos \theta) e^{i m \phi}
where:
Copy link
Member

@f-brinkmann f-brinkmann Mar 27, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The math below is not displayed correctly in the docs. Probably you need to use the :math: tag in Sphinx

Comment on lines +265 to +267
- $P_{nm}$ is the associated Legendre function
- $N_{nm}$ is the normalization term
- $CS_m$ is the Condon-Shortley phase term
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we link somewhere for an analytical description of these? For the last, we could also write CS_m = (-1)^m directly here.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copy link
Member

@f-brinkmann f-brinkmann left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Found only a few minor things. I generally wonder where we give a complete analytical description of the SHs? Is that given in the pull related to the SH-class?

Copy link
Member

@ahms5 ahms5 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you for taking care! just checked the docs

@@ -159,6 +159,8 @@ def _spherical_hankel_derivative(n, z, kind):

def spherical_harmonic(n, m, theta, phi):
"""The spherical harmonics of order n and degree m.
The spherical harmonic functions are fully normalized (N3D) and
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
The spherical harmonic functions are fully normalized (N3D) and
The spherical harmonic functions are fully normalized (N3D) and

very picky suggestion^^

Comment on lines +205 to +206
the AmbiX phase convention which does not include the Condon-Shortley
phase [#]_.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
the AmbiX phase convention which does not include the Condon-Shortley
phase [#]_.
the AmbiX phase convention, which does not include the Condon-Shortley
phase [#]_.


The spherical harmonic functions are fully normalized (N3D) and include the
Condon-Shotley phase term :math:`(-1)^m` [#]_, [#]_.
See also :func:`spherical_harmonic_basis_real`.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
See also :func:`spherical_harmonic_basis_real`.
See [#]_ and [#]_ for details.
See also :py:func:`spherical_harmonic_basis_real`.

we need to refer the references, otherwise the doc build on ci will fail

Y_n^m(\theta, \phi) = \sqrt{\frac{2n+1}{4\pi}
\frac{(n-m)!}{(n+m)!}} P_n^m(\cos \theta) e^{i m \phi}
where:
- $n$ is the degree
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
- $n$ is the degree
- :math:`n` is the degree

like this, the empty line should also fix the bulletpoints

(FuMa is only supported up to 3rd order)
phase_convention : string or None, optional
Whether to include the Condon-Shortley phase term.
The default is 'Condon-Shortley'.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

we should test if no other string than 'Condon-Shortley' is used


The spherical harmonic functions are fully normalized (N3D) and follow
the AmbiX phase convention [#]_.
See also :func:`spherical_harmonic_basis`.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
See also :func:`spherical_harmonic_basis`.
See also :py:func:`spherical_harmonic_basis`.

(FuMa is only supported up to 3rd order)
phase_convention : string or None, optional
Whether to include the Condon-Shortley phase term.
The default is None.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

does it make sence to have different defaults for the "same" parameter?
we should also highlight the possible opitons.
if we wouldnt add more options in future we could also think about a boolean input ike apply_cs_phase

@@ -471,9 +567,21 @@ def spherical_harmonic_basis_gradient_real(n_max, coords):
----------
n_max : int
Spherical harmonic order
coordinates : :doc:`pf.Coordinates <pyfar:classes/pyfar.coordinates>`
Coordinate object with sampling points for which the basis matrix is
coordinates : :doc:`pf.Coordinates <pyfar:classes/pyfar.coordinates>` or
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

see above

def spherical_harmonic_basis_gradient_real(n_max, coordinates,
normalization="n3d",
channel_convention="acn",
phase_convention=None):
r"""
Calulcates the unit sphere gradients of the real valued spherical hamonics.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Calculates harmonics*

(FuMa is only supported up to 3rd order)
phase_convention : string or None, optional
Whether to include the Condon-Shortley phase term.
The default is None.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

see above

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
Status: Require review
Development

Successfully merging this pull request may close these issues.

5 participants