From 068c31949968e2cac0f9d445ded95dade1320a84 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Tue, 31 Mar 2026 20:27:59 +0000
Subject: [PATCH 1/6] Add SphericalHarmonicsFilter with
 convolve/multiply/rotate/from_distribution_numerical_fast

Agent-Logs-Url: https://github.com/FlorianPfaff/PyRecEst/sessions/a02e218a-d9bc-4d03-9c38-8776e308de70

Co-authored-by: FlorianPfaff <6773539+FlorianPfaff@users.noreply.github.com>
---
 ...pherical_harmonics_distribution_complex.py | 260 ++++++++++++++++++
 pyrecest/filters/__init__.py                  |   2 +
 .../filters/spherical_harmonics_filter.py     | 229 +++++++++++++++
 pyrecest/filters/von_mises_fisher_filter.py   |   8 +
 .../test_spherical_harmonics_filter.py        | 112 ++++++++
 5 files changed, 611 insertions(+)
 create mode 100644 pyrecest/filters/spherical_harmonics_filter.py
 create mode 100644 pyrecest/tests/filters/test_spherical_harmonics_filter.py

diff --git a/pyrecest/distributions/hypersphere_subset/spherical_harmonics_distribution_complex.py b/pyrecest/distributions/hypersphere_subset/spherical_harmonics_distribution_complex.py
index 03ed421a7..141e9cecd 100644
--- a/pyrecest/distributions/hypersphere_subset/spherical_harmonics_distribution_complex.py
+++ b/pyrecest/distributions/hypersphere_subset/spherical_harmonics_distribution_complex.py
@@ -1,3 +1,6 @@
+import warnings
+
+import numpy as np
 import scipy
 
 # pylint: disable=redefined-builtin,no-name-in-module,no-member
@@ -190,3 +193,260 @@ def imag_part(phi, theta, n, m):
                 coeff_mat[n, m + n] = real_integral + 1j * imag_integral
 
         return SphericalHarmonicsDistributionComplex(coeff_mat, transformation)
+
+    # ------------------------------------------------------------------
+    # pyshtools-based helper methods
+    # ------------------------------------------------------------------
+
+    @staticmethod
+    def _coeff_mat_to_pysh(coeff_mat, degree):
+        """Convert our coeff_mat to a pyshtools SHComplexCoeffs object."""
+        import pyshtools as pysh  # pylint: disable=import-error
+
+        clm = pysh.SHCoeffs.from_zeros(
+            degree, kind="complex", normalization="ortho", csphase=-1
+        )
+        for n in range(degree + 1):
+            for m in range(0, n + 1):
+                clm.coeffs[0, n, m] = coeff_mat[n, n + m]
+            for m in range(1, n + 1):
+                clm.coeffs[1, n, m] = coeff_mat[n, n - m]
+        return clm
+
+    @staticmethod
+    def _pysh_to_coeff_mat(clm, degree):
+        """Convert a pyshtools SHComplexCoeffs object to our coeff_mat."""
+        coeff_mat = np.zeros((degree + 1, 2 * degree + 1), dtype=complex)
+        max_n = min(clm.lmax, degree)
+        for n in range(max_n + 1):
+            for m in range(0, n + 1):
+                coeff_mat[n, n + m] = clm.coeffs[0, n, m]
+            for m in range(1, n + 1):
+                coeff_mat[n, n - m] = clm.coeffs[1, n, m]
+        return coeff_mat
+
+    @staticmethod
+    def _get_dh_grid_cartesian(degree):
+        """Return (x, y, z) flat arrays and grid_shape for the DH grid at *degree*."""
+        import pyshtools as pysh  # pylint: disable=import-error
+
+        dummy = pysh.SHCoeffs.from_zeros(
+            degree, kind="complex", normalization="ortho", csphase=-1
+        )
+        grid = dummy.expand(grid="DH", extend=False)
+        lats, lons = grid.lats(), grid.lons()
+        lon_mesh, lat_mesh = np.meshgrid(lons, lats)
+        theta = np.radians(90.0 - lat_mesh)  # colatitude in radians
+        phi = np.radians(lon_mesh)  # azimuth in radians
+        x_c = np.sin(theta) * np.cos(phi)
+        y_c = np.sin(theta) * np.sin(phi)
+        z_c = np.cos(theta)
+        return x_c.ravel(), y_c.ravel(), z_c.ravel(), theta.shape
+
+    def _eval_on_grid(self, target_degree=None):
+        """Evaluate this SHD on the DH grid at *target_degree* (defaults to own degree).
+
+        The DH grid is expanded at *target_degree* so that higher-frequency grids
+        can be used for intermediate computations (e.g. when squaring a degree-L
+        function that has degree 2L).
+        Returns a real 2-D numpy array of shape (nlat, nlon).
+        """
+        import pyshtools as pysh  # pylint: disable=import-error
+
+        degree = self.coeff_mat.shape[0] - 1
+        if target_degree is None:
+            target_degree = degree
+
+        # Pad our coefficients into a pysh object at target_degree
+        clm_full = pysh.SHCoeffs.from_zeros(
+            target_degree, kind="complex", normalization="ortho", csphase=-1
+        )
+        min_deg = min(degree, target_degree)
+        clm_own = self._coeff_mat_to_pysh(self.coeff_mat, min_deg)
+        clm_full.coeffs[0, : min_deg + 1, : min_deg + 1] = clm_own.coeffs[
+            0, : min_deg + 1, : min_deg + 1
+        ]
+        clm_full.coeffs[1, : min_deg + 1, : min_deg + 1] = clm_own.coeffs[
+            1, : min_deg + 1, : min_deg + 1
+        ]
+        grid = clm_full.expand(grid="DH", extend=False)
+        return grid.data.real
+
+    @staticmethod
+    def _fit_from_grid(grid_vals_real, degree, transformation):
+        """Fit SH coefficients to real-valued grid values on a DH grid.
+
+        *grid_vals_real* is a 2D numpy array with shape matching the DH grid for
+        *degree*.  Returns a new :class:`SphericalHarmonicsDistributionComplex`.
+        """
+        import pyshtools as pysh  # pylint: disable=import-error
+
+        grid_obj = pysh.SHGrid.from_array(
+            grid_vals_real.astype(complex), grid="DH"
+        )
+        clm = grid_obj.expand(lmax_calc=degree, normalization="ortho", csphase=-1)
+        coeff_mat = SphericalHarmonicsDistributionComplex._pysh_to_coeff_mat(clm, degree)
+        shd = SphericalHarmonicsDistributionComplex(coeff_mat, transformation)
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            shd.normalize_in_place()
+        return shd
+
+    # ------------------------------------------------------------------
+    # Public numerical methods
+    # ------------------------------------------------------------------
+
+    @staticmethod
+    def from_distribution_numerical_fast(dist, degree, transformation="identity"):
+        """Approximate *dist* by a degree-*degree* SHD using a DH grid.
+
+        This is faster than :meth:`from_distribution_via_integral` because it
+        uses the discrete spherical-harmonic transform instead of numerical
+        quadrature.
+        """
+        if transformation not in ("identity", "sqrt"):
+            raise ValueError(f"Unsupported transformation: '{transformation}'")
+
+        x_c, y_c, z_c, grid_shape = (
+            SphericalHarmonicsDistributionComplex._get_dh_grid_cartesian(degree)
+        )
+        xs = np.column_stack([x_c, y_c, z_c])
+        fvals = np.asarray(dist.pdf(xs), dtype=float).reshape(grid_shape)
+
+        if transformation == "sqrt":
+            fvals = np.sqrt(np.maximum(fvals, 0.0))
+
+        return SphericalHarmonicsDistributionComplex._fit_from_grid(
+            fvals, degree, transformation
+        )
+
+    def convolve(self, other):
+        """Spherical convolution with a *zonal* distribution *other*.
+
+        For the ``'identity'`` transformation the standard frequency-domain
+        formula is used (exact for bandlimited functions).  For the ``'sqrt'``
+        transformation a grid-based approach with a 2× finer intermediate grid
+        is used so that squaring the sqrt functions introduces no aliasing.
+        """
+        assert isinstance(
+            other, SphericalHarmonicsDistributionComplex
+        ), "other must be a SphericalHarmonicsDistributionComplex"
+
+        degree = self.coeff_mat.shape[0] - 1
+
+        if self.transformation == "identity" and other.transformation == "identity":
+            # Direct frequency-domain formula: h_{l,m} = sqrt(4π/(2l+1)) * f_{l,m} * g_{l,0}
+            h_lm = np.zeros_like(np.asarray(self.coeff_mat))
+            for l in range(degree + 1):
+                factor = (
+                    np.sqrt(4.0 * np.pi / (2 * l + 1))
+                    * np.asarray(other.coeff_mat[l, l])
+                )
+                for m in range(-l, l + 1):
+                    h_lm[l, l + m] = (
+                        factor * np.asarray(self.coeff_mat[l, l + m])
+                    )
+            result = SphericalHarmonicsDistributionComplex(h_lm, "identity")
+            with warnings.catch_warnings():
+                warnings.simplefilter("ignore")
+                result.normalize_in_place()
+            return result
+
+        if self.transformation == "sqrt" and other.transformation == "sqrt":
+            # Use a grid twice as fine to avoid aliasing when squaring.
+            degree_fine = 2 * degree
+
+            # Recover p = f^2 and q = g^2 on the fine grid, then SHT to degree.
+            f_grid = self._eval_on_grid(target_degree=degree_fine)
+            p_grid = f_grid**2
+
+            g_grid = other._eval_on_grid(target_degree=degree_fine)  # pylint: disable=protected-access
+            q_grid = g_grid**2
+
+            import pyshtools as pysh  # pylint: disable=import-error
+
+            def _grid_to_coeff(grid_vals):
+                grid_obj = pysh.SHGrid.from_array(
+                    grid_vals.astype(complex), grid="DH"
+                )
+                clm = grid_obj.expand(
+                    lmax_calc=degree, normalization="ortho", csphase=-1
+                )
+                return self._pysh_to_coeff_mat(clm, degree)
+
+            p_lm = _grid_to_coeff(p_grid)
+            q_lm = _grid_to_coeff(q_grid)
+
+            # Convolution formula on the identity coefficients
+            r_lm = np.zeros_like(p_lm)
+            for l in range(degree + 1):
+                factor = np.sqrt(4.0 * np.pi / (2 * l + 1)) * q_lm[l, l]
+                for m in range(-l, l + 1):
+                    r_lm[l, l + m] = factor * p_lm[l, l + m]
+
+            # Evaluate r on the standard DH grid, take sqrt, refit
+            r_shd_id = SphericalHarmonicsDistributionComplex(r_lm, "identity")
+            r_grid = r_shd_id._eval_on_grid()
+            sqrt_r_grid = np.sqrt(np.maximum(r_grid, 0.0))
+
+            return SphericalHarmonicsDistributionComplex._fit_from_grid(
+                sqrt_r_grid, degree, "sqrt"
+            )
+
+        raise ValueError(
+            "convolve: mixed transformations are not supported. "
+            "Both self and other must use the same transformation."
+        )
+
+    def multiply(self, other):
+        """Pointwise multiplication in physical space (Bayesian update step).
+
+        Works for both ``'identity'`` and ``'sqrt'`` transformations.
+        For ``'sqrt'``: ``sqrt(p) * sqrt(q) = sqrt(p*q)`` which is the correct
+        sqrt-transformed product density.
+        """
+        assert isinstance(
+            other, SphericalHarmonicsDistributionComplex
+        ), "other must be a SphericalHarmonicsDistributionComplex"
+        assert self.transformation == other.transformation, (
+            "multiply: both distributions must use the same transformation"
+        )
+
+        degree = self.coeff_mat.shape[0] - 1
+
+        f_grid = self._eval_on_grid()
+        g_grid = other._eval_on_grid()
+
+        h_grid = f_grid * g_grid
+
+        return SphericalHarmonicsDistributionComplex._fit_from_grid(
+            h_grid, degree, self.transformation
+        )
+
+    def rotate(self, alpha, beta, gamma):
+        """Rotate the distribution by ZYZ Euler angles (in radians).
+
+        Parameters
+        ----------
+        alpha : float
+            First rotation angle around Z (radians).
+        beta : float
+            Second rotation angle around Y (radians).
+        gamma : float
+            Third rotation angle around Z (radians).
+        """
+        import pyshtools as pysh  # pylint: disable=import-error
+
+        degree = self.coeff_mat.shape[0] - 1
+        clm = self._coeff_mat_to_pysh(self.coeff_mat, degree)
+        clm_rot = clm.rotate(
+            np.degrees(alpha),
+            np.degrees(beta),
+            np.degrees(gamma),
+            degrees=True,
+            body=True,
+        )
+        coeff_mat_rot = self._pysh_to_coeff_mat(clm_rot, degree)
+        return SphericalHarmonicsDistributionComplex(
+            coeff_mat_rot, self.transformation
+        )
diff --git a/pyrecest/filters/__init__.py b/pyrecest/filters/__init__.py
index 914306392..cbe26252a 100644
--- a/pyrecest/filters/__init__.py
+++ b/pyrecest/filters/__init__.py
@@ -4,6 +4,7 @@
 from .hypertoroidal_particle_filter import HypertoroidalParticleFilter
 from .kalman_filter import KalmanFilter
 from .manifold_mixins import EuclideanFilterMixin, HypertoroidalFilterMixin
+from .von_mises_fisher_filter import VonMisesFisherFilter
 
 __all__ = [
     "AbstractFilter",
@@ -13,4 +14,5 @@
     "HypertoroidalParticleFilter",
     "KalmanFilter",
     "EuclideanParticleFilter",
+    "VonMisesFisherFilter",
 ]
diff --git a/pyrecest/filters/spherical_harmonics_filter.py b/pyrecest/filters/spherical_harmonics_filter.py
new file mode 100644
index 000000000..f24f42bac
--- /dev/null
+++ b/pyrecest/filters/spherical_harmonics_filter.py
@@ -0,0 +1,229 @@
+import copy
+import warnings
+
+import numpy as np
+
+from pyrecest.distributions.hypersphere_subset.spherical_harmonics_distribution_complex import (
+    SphericalHarmonicsDistributionComplex,
+)
+
+from .abstract_filter import AbstractFilter
+from .manifold_mixins import HypersphericalFilterMixin
+
+
+class SphericalHarmonicsFilter(AbstractFilter, HypersphericalFilterMixin):
+    """Filter on the unit sphere using spherical harmonic representations.
+
+    Supports both the ``'identity'`` transformation (coefficients represent
+    the density directly) and the ``'sqrt'`` transformation (coefficients
+    represent the square-root of the density).
+
+    References
+    ----------
+    Florian Pfaff, Gerhard Kurz, and Uwe D. Hanebeck,
+    "Filtering on the Unit Sphere Using Spherical Harmonics",
+    Proceedings of the 2017 IEEE International Conference on Multisensor
+    Fusion and Integration for Intelligent Systems (MFI 2017),
+    Daegu, Korea, November 2017.
+    """
+
+    def __init__(self, degree, transformation="identity"):
+        HypersphericalFilterMixin.__init__(self)
+        coeff_mat = np.zeros((degree + 1, 2 * degree + 1), dtype=complex)
+        if transformation == "identity":
+            coeff_mat[0, 0] = 1.0 / np.sqrt(4.0 * np.pi)
+        elif transformation == "sqrt":
+            coeff_mat[0, 0] = 1.0
+        else:
+            raise ValueError(f"Unknown transformation: '{transformation}'")
+        initial_state = SphericalHarmonicsDistributionComplex(
+            coeff_mat, transformation
+        )
+        AbstractFilter.__init__(self, initial_state)
+
+    # ------------------------------------------------------------------
+    # filter_state / state property
+    # ------------------------------------------------------------------
+
+    @property
+    def filter_state(self):
+        return self._filter_state
+
+    @filter_state.setter
+    def filter_state(self, new_state):
+        assert isinstance(
+            new_state, SphericalHarmonicsDistributionComplex
+        ), "filter_state must be a SphericalHarmonicsDistributionComplex"
+        self._filter_state = copy.deepcopy(new_state)
+
+    @property
+    def state(self):
+        """Alias for :attr:`filter_state`."""
+        return self._filter_state
+
+    # ------------------------------------------------------------------
+    # Public interface methods
+    # ------------------------------------------------------------------
+
+    def set_state(self, state):
+        """Set the filter state with optional warnings about mismatches."""
+        assert isinstance(
+            state, SphericalHarmonicsDistributionComplex
+        ), "state must be a SphericalHarmonicsDistributionComplex"
+        if self._filter_state.transformation != state.transformation:
+            warnings.warn(
+                "setState:transDiffer: New density is transformed differently.",
+                stacklevel=2,
+            )
+        if self._filter_state.coeff_mat.shape != state.coeff_mat.shape:
+            warnings.warn(
+                "setState:noOfCoeffsDiffer: New density has different number of "
+                "coefficients.",
+                stacklevel=2,
+            )
+        self._filter_state = copy.deepcopy(state)
+
+    def get_estimate(self):
+        """Return the current filter state."""
+        return self._filter_state
+
+    def get_estimate_mean(self):
+        """Return the mean direction of the current filter state."""
+        return self._filter_state.mean_direction()
+
+    # ------------------------------------------------------------------
+    # Prediction
+    # ------------------------------------------------------------------
+
+    def predict_identity(self, sys_noise):
+        """Predict via spherical convolution with a *zonal* system noise SHD.
+
+        Parameters
+        ----------
+        sys_noise : SphericalHarmonicsDistributionComplex
+            Must be a zonal distribution (rotationally symmetric around the
+            z-axis) in the same transformation as the filter state.
+        """
+        assert isinstance(
+            sys_noise, SphericalHarmonicsDistributionComplex
+        ), "sys_noise must be a SphericalHarmonicsDistributionComplex"
+        if (
+            self._filter_state.transformation == "sqrt"
+            and sys_noise.transformation == "identity"
+        ):
+            state_degree = self._filter_state.coeff_mat.shape[0] - 1
+            noise_degree = sys_noise.coeff_mat.shape[0] - 1
+            assert 2 * state_degree == noise_degree, (
+                "If the sqrt variant is used and sys_noise is given in "
+                "identity form, sys_noise should have degree 2 * state_degree."
+            )
+        self._filter_state = self._filter_state.convolve(sys_noise)
+
+    # ------------------------------------------------------------------
+    # Update
+    # ------------------------------------------------------------------
+
+    def update_identity(self, meas_noise, z):
+        """Update by multiplying the state with a (possibly rotated) noise SHD.
+
+        *meas_noise* should be a zonal SHD with its axis of symmetry along
+        [0, 0, 1].  If the measurement *z* differs from [0, 0, 1] the noise is
+        rotated to align with *z* before the multiplication.
+
+        Parameters
+        ----------
+        meas_noise : SphericalHarmonicsDistributionComplex
+            Zonal measurement noise (axis along [0, 0, 1]).
+        z : array-like, shape (3,)
+            Measurement direction on the unit sphere.
+        """
+        assert isinstance(
+            meas_noise, SphericalHarmonicsDistributionComplex
+        ), "meas_noise must be a SphericalHarmonicsDistributionComplex"
+        z = np.asarray(z, dtype=float).ravel()
+        z_norm = np.linalg.norm(z)
+        if z_norm > 1e-6 and np.linalg.norm(z - np.array([0.0, 0.0, 1.0])) > 1e-6:
+            warnings.warn(
+                "SphericalHarmonicsFilter:rotationRequired: "
+                "Performance may be low for z != [0, 0, 1]. "
+                "Using update_nonlinear may yield faster results.",
+                stacklevel=2,
+            )
+            phi = np.arctan2(z[1], z[0])  # azimuth
+            theta = np.arccos(
+                np.clip(z[2] / z_norm, -1.0, 1.0)
+            )  # colatitude
+            meas_noise = meas_noise.rotate(0.0, theta, phi)
+        self._filter_state = self._filter_state.multiply(meas_noise)
+
+    def update_nonlinear(self, likelihood, z):
+        """Nonlinear Bayesian update via a likelihood function.
+
+        Parameters
+        ----------
+        likelihood : callable
+            ``likelihood(z, pts)`` where *pts* is a ``(3, N)`` matrix of
+            Cartesian coordinates on the unit sphere and the return value is a
+            length-N array of likelihood values.
+        z : array-like
+            Measurement (passed through to *likelihood* unchanged).
+        """
+        self._update_nonlinear_impl([likelihood], [z])
+
+    def update_nonlinear_multiple(self, likelihoods, measurements):
+        """Nonlinear update using a list of likelihood functions simultaneously.
+
+        Parameters
+        ----------
+        likelihoods : list of callables
+            Each element is a likelihood function as described in
+            :meth:`update_nonlinear`.
+        measurements : list of array-like
+            Corresponding measurements.
+        """
+        assert len(likelihoods) == len(
+            measurements
+        ), "likelihoods and measurements must have the same length"
+        self._update_nonlinear_impl(likelihoods, measurements)
+
+    # ------------------------------------------------------------------
+    # Internal helpers
+    # ------------------------------------------------------------------
+
+    def _update_nonlinear_impl(self, likelihoods, measurements):
+        """Shared implementation for single and multiple nonlinear updates."""
+        import pyshtools as pysh  # pylint: disable=import-error
+
+        degree = self._filter_state.coeff_mat.shape[0] - 1
+
+        # DH grid coordinates
+        x_c, y_c, z_c, grid_shape = (
+            SphericalHarmonicsDistributionComplex._get_dh_grid_cartesian(degree)
+        )
+        # (3, N) matrix for likelihood calls
+        grid_pts = np.stack([x_c, y_c, z_c], axis=0)
+
+        # Evaluate current state on the DH grid
+        fval_curr = self._filter_state._eval_on_grid()  # pylint: disable=protected-access
+
+        # Accumulate likelihood values over all (likelihood, measurement) pairs
+        likelihood_vals = np.ones(grid_shape, dtype=float)
+        for lk, zk in zip(likelihoods, measurements):
+            lv = np.asarray(lk(zk, grid_pts), dtype=float).reshape(grid_shape)
+            likelihood_vals *= lv
+
+        # Scale factor (improves numerical conditioning near zero)
+        scale = float(2 ** len(likelihoods))
+
+        if self._filter_state.transformation == "identity":
+            fval_new = scale * fval_curr * likelihood_vals
+        elif self._filter_state.transformation == "sqrt":
+            fval_new = scale * fval_curr * np.sqrt(np.maximum(likelihood_vals, 0.0))
+        else:
+            raise ValueError(
+                f"Unsupported transformation: '{self._filter_state.transformation}'"
+            )
+
+        self._filter_state = SphericalHarmonicsDistributionComplex._fit_from_grid(
+            fval_new, degree, self._filter_state.transformation
+        )
diff --git a/pyrecest/filters/von_mises_fisher_filter.py b/pyrecest/filters/von_mises_fisher_filter.py
index 2ea742143..6e274da9a 100644
--- a/pyrecest/filters/von_mises_fisher_filter.py
+++ b/pyrecest/filters/von_mises_fisher_filter.py
@@ -24,6 +24,14 @@ def filter_state(self, filter_state):
         ), "filter_state must be an instance of VonMisesFisherDistribution."
         self._filter_state = filter_state
 
+    def set_state(self, state):
+        """Set the filter state."""
+        self.filter_state = state
+
+    def get_estimate_mean(self):
+        """Return the mean direction of the current filter state."""
+        return self.filter_state.mean_direction()
+
     def predict_identity(self, sys_noise):
         """
         State prediction via mulitiplication. Provide zonal density for update
diff --git a/pyrecest/tests/filters/test_spherical_harmonics_filter.py b/pyrecest/tests/filters/test_spherical_harmonics_filter.py
new file mode 100644
index 000000000..43e9c15a9
--- /dev/null
+++ b/pyrecest/tests/filters/test_spherical_harmonics_filter.py
@@ -0,0 +1,112 @@
+import unittest
+import numpy as np
+from scipy.stats import norm
+
+from pyrecest.distributions import VonMisesFisherDistribution, SphericalHarmonicsDistributionComplex
+from pyrecest.filters import VonMisesFisherFilter
+from pyrecest.filters.spherical_harmonics_filter import SphericalHarmonicsFilter
+
+class SphericalHarmonicsFilterTest(unittest.TestCase):
+
+    def test_update_identity(self):
+        for transformation in ['identity', 'sqrt']:
+            shd_filter = SphericalHarmonicsFilter(30, transformation)
+            vmf_filter = VonMisesFisherFilter()
+
+            vmf1 = VonMisesFisherDistribution(np.array([0, 1, 0]), 1)
+            vmf2 = VonMisesFisherDistribution(np.array([0, 0, 1]), 0.1)
+            
+            shd1 = SphericalHarmonicsDistributionComplex.from_distribution_numerical_fast(vmf1, 30, transformation)
+            shd2 = SphericalHarmonicsDistributionComplex.from_distribution_numerical_fast(vmf2, 30, transformation)
+
+            vmf_filter.set_state(vmf1)
+            vmf_filter.update_identity(vmf2, np.array([1, 0, 0]))
+
+            shd_filter.set_state(shd1)
+            shd_filter.update_identity(shd2, np.array([1, 0, 0]))
+
+            self.assertTrue(np.allclose(vmf_filter.get_estimate_mean(), shd_filter.get_estimate_mean(), atol=1e-10))
+
+    def test_update_using_likelihood(self):
+        np.random.seed(1)
+
+        pos_true = -1 / np.sqrt(3) * np.ones(3)
+
+        # Generate measurements according to truncated gaussian along x, y, and z axis
+        sigma_x = 0.3
+        sigma_y = 0.3
+        sigma_z = 0.3
+        
+        meas_x = self._generate_truncated_normals(pos_true[0], sigma_x, 5)
+        meas_y = self._generate_truncated_normals(pos_true[1], sigma_y, 5)
+        meas_z = self._generate_truncated_normals(pos_true[2], sigma_z, 5)
+
+        for transformation in ['identity', 'sqrt']:
+            sh_filter = SphericalHarmonicsFilter(11, transformation)
+
+            for x in meas_x:
+                sh_filter.update_nonlinear(lambda z, x: norm.pdf(z[0], x[0], sigma_x), np.array([x, 0, 0]))
+            for y in meas_y:
+                sh_filter.update_nonlinear(lambda z, x: norm.pdf(z[1], x[1], sigma_y), np.array([0, y, 0]))
+            for z in meas_z:
+                sh_filter.update_nonlinear(lambda z, x: norm.pdf(z[2], x[2], sigma_z), np.array([0, 0, z]))
+            
+            self.assertTrue(np.allclose(sh_filter.get_estimate_mean(), pos_true, atol=0.3))
+
+    def test_update_using_likelihood_multiple(self):
+        for transformation in ['identity', 'sqrt']:
+            sh_filter1 = SphericalHarmonicsFilter(10, transformation)
+            sh_filter2 = SphericalHarmonicsFilter(10, transformation)
+
+            sigma_x = 0.3
+            sigma_y = 0.3
+            sigma_z = 0.3
+            
+            sh_filter1.update_nonlinear(lambda z, x: norm.pdf(z[0], x[0], sigma_x), np.array([-1 / np.sqrt(3), 0, 0]))
+            sh_filter1.update_nonlinear(lambda z, x: norm.pdf(z[1], x[1], sigma_y), np.array([0, -1 / np.sqrt(3), 0]))
+            sh_filter1.update_nonlinear(lambda z, x: norm.pdf(z[2], x[2], sigma_z), np.array([0, 0, -1 / np.sqrt(3)]))
+
+            sh_filter2.update_nonlinear_multiple(
+                [
+                    lambda z, x: norm.pdf(z[0], x[0], sigma_x),
+                    lambda z, x: norm.pdf(z[1], x[1], sigma_y),
+                    lambda z, x: norm.pdf(z[2], x[2], sigma_z)
+                ],
+                [
+                    np.array([-1 / np.sqrt(3), 0, 0]),
+                    np.array([0, -1 / np.sqrt(3), 0]),
+                    np.array([0, 0, -1 / np.sqrt(3)])
+                ]
+            )
+
+            self.assertTrue(np.allclose(sh_filter2.get_estimate_mean(), sh_filter1.get_estimate_mean(), atol=1e-5))
+
+    def test_prediction_sqrt_vs_id(self):
+        degree = 21
+        density_init = VonMisesFisherDistribution(np.array([1, 1, 0]) / np.sqrt(2), 2)
+        sys_noise = VonMisesFisherDistribution(np.array([0, 0, 1]), 1)
+
+        shd_init_id = SphericalHarmonicsDistributionComplex.from_distribution_numerical_fast(density_init, degree, 'identity')
+        shd_init_sqrt = SphericalHarmonicsDistributionComplex.from_distribution_numerical_fast(density_init, degree, 'sqrt')
+        shd_noise_id = SphericalHarmonicsDistributionComplex.from_distribution_numerical_fast(sys_noise, degree, 'identity')
+        shd_noise_sqrt = SphericalHarmonicsDistributionComplex.from_distribution_numerical_fast(sys_noise, degree, 'sqrt')
+
+        sh_filter1 = SphericalHarmonicsFilter(degree, 'identity')
+        sh_filter2 = SphericalHarmonicsFilter(degree, 'sqrt')
+
+        sh_filter1.set_state(shd_init_id)
+        sh_filter2.set_state(shd_init_sqrt)
+
+        sh_filter1.predict_identity(shd_noise_id)
+        sh_filter2.predict_identity(shd_noise_sqrt)
+
+        np.testing.assert_allclose(sh_filter1.state.total_variation_distance_numerical(sh_filter2.state), 0, atol=5e-15)
+
+    # Helper function to generate truncated normals
+    def _generate_truncated_normals(self, mu, sigma, n):
+        samples = []
+        while len(samples) < n:
+            sample = np.random.normal(mu, sigma)
+            if -1 <= sample <= 1:
+                samples.append(sample)
+        return samples

From d4c3a4d304c2e9123645184b74267f8145a5ca50 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Tue, 31 Mar 2026 20:31:14 +0000
Subject: [PATCH 2/6] Address code review: remove unnecessary np.asarray, fix
 z-check, add scale factor comment

Agent-Logs-Url: https://github.com/FlorianPfaff/PyRecEst/sessions/a02e218a-d9bc-4d03-9c38-8776e308de70

Co-authored-by: FlorianPfaff <6773539+FlorianPfaff@users.noreply.github.com>
---
 .../spherical_harmonics_distribution_complex.py        | 10 ++++------
 pyrecest/filters/spherical_harmonics_filter.py         | 10 ++++++++--
 2 files changed, 12 insertions(+), 8 deletions(-)

diff --git a/pyrecest/distributions/hypersphere_subset/spherical_harmonics_distribution_complex.py b/pyrecest/distributions/hypersphere_subset/spherical_harmonics_distribution_complex.py
index 141e9cecd..4342ecfeb 100644
--- a/pyrecest/distributions/hypersphere_subset/spherical_harmonics_distribution_complex.py
+++ b/pyrecest/distributions/hypersphere_subset/spherical_harmonics_distribution_complex.py
@@ -336,16 +336,14 @@ def convolve(self, other):
 
         if self.transformation == "identity" and other.transformation == "identity":
             # Direct frequency-domain formula: h_{l,m} = sqrt(4π/(2l+1)) * f_{l,m} * g_{l,0}
-            h_lm = np.zeros_like(np.asarray(self.coeff_mat))
+            h_lm = np.zeros_like(self.coeff_mat)
             for l in range(degree + 1):
                 factor = (
                     np.sqrt(4.0 * np.pi / (2 * l + 1))
-                    * np.asarray(other.coeff_mat[l, l])
+                    * other.coeff_mat[l, l]
                 )
                 for m in range(-l, l + 1):
-                    h_lm[l, l + m] = (
-                        factor * np.asarray(self.coeff_mat[l, l + m])
-                    )
+                    h_lm[l, l + m] = factor * self.coeff_mat[l, l + m]
             result = SphericalHarmonicsDistributionComplex(h_lm, "identity")
             with warnings.catch_warnings():
                 warnings.simplefilter("ignore")
@@ -360,7 +358,7 @@ def convolve(self, other):
             f_grid = self._eval_on_grid(target_degree=degree_fine)
             p_grid = f_grid**2
 
-            g_grid = other._eval_on_grid(target_degree=degree_fine)  # pylint: disable=protected-access
+            g_grid = other._eval_on_grid(target_degree=degree_fine)
             q_grid = g_grid**2
 
             import pyshtools as pysh  # pylint: disable=import-error
diff --git a/pyrecest/filters/spherical_harmonics_filter.py b/pyrecest/filters/spherical_harmonics_filter.py
index f24f42bac..6f0ba7615 100644
--- a/pyrecest/filters/spherical_harmonics_filter.py
+++ b/pyrecest/filters/spherical_harmonics_filter.py
@@ -142,7 +142,10 @@ def update_identity(self, meas_noise, z):
         ), "meas_noise must be a SphericalHarmonicsDistributionComplex"
         z = np.asarray(z, dtype=float).ravel()
         z_norm = np.linalg.norm(z)
-        if z_norm > 1e-6 and np.linalg.norm(z - np.array([0.0, 0.0, 1.0])) > 1e-6:
+        not_near_north_pole = (
+            np.abs(z[0]) > 1e-6 or np.abs(z[1]) > 1e-6 or np.abs(z[2] - 1.0) > 1e-6
+        )
+        if z_norm > 1e-6 and not_near_north_pole:
             warnings.warn(
                 "SphericalHarmonicsFilter:rotationRequired: "
                 "Performance may be low for z != [0, 0, 1]. "
@@ -212,7 +215,10 @@ def _update_nonlinear_impl(self, likelihoods, measurements):
             lv = np.asarray(lk(zk, grid_pts), dtype=float).reshape(grid_shape)
             likelihood_vals *= lv
 
-        # Scale factor (improves numerical conditioning near zero)
+        # Scale factor: multiplying by 2^n keeps values away from zero so that
+        # the SHT fit and subsequent normalisation remain numerically stable
+        # when the product likelihood is very small (e.g. many weak likelihoods).
+        # This factor is divided out implicitly by the normalisation step.
         scale = float(2 ** len(likelihoods))
 
         if self._filter_state.transformation == "identity":

From c527e6df6cde3b90a527bc1c04c99810e54fbcad Mon Sep 17 00:00:00 2001
From: Florian Pfaff <pfaff@ias.uni-stuttgart.de>
Date: Wed, 1 Apr 2026 01:52:18 +0200
Subject: [PATCH 3/6] Support all backends for spherical harmonics filter

---
 .../filters/spherical_harmonics_filter.py     | 26 +++++++++----------
 1 file changed, 13 insertions(+), 13 deletions(-)

diff --git a/pyrecest/filters/spherical_harmonics_filter.py b/pyrecest/filters/spherical_harmonics_filter.py
index 6f0ba7615..2df579aac 100644
--- a/pyrecest/filters/spherical_harmonics_filter.py
+++ b/pyrecest/filters/spherical_harmonics_filter.py
@@ -1,7 +1,7 @@
 import copy
 import warnings
 
-import numpy as np
+from pyrecest.backend import zeros, sqrt, pi, linalg, abs, stack, arccos, arctan2, clip, array, ones
 
 from pyrecest.distributions.hypersphere_subset.spherical_harmonics_distribution_complex import (
     SphericalHarmonicsDistributionComplex,
@@ -29,9 +29,9 @@ class SphericalHarmonicsFilter(AbstractFilter, HypersphericalFilterMixin):
 
     def __init__(self, degree, transformation="identity"):
         HypersphericalFilterMixin.__init__(self)
-        coeff_mat = np.zeros((degree + 1, 2 * degree + 1), dtype=complex)
+        coeff_mat = zeros((degree + 1, 2 * degree + 1), dtype=complex)
         if transformation == "identity":
-            coeff_mat[0, 0] = 1.0 / np.sqrt(4.0 * np.pi)
+            coeff_mat[0, 0] = 1.0 / sqrt(4.0 * pi)
         elif transformation == "sqrt":
             coeff_mat[0, 0] = 1.0
         else:
@@ -140,10 +140,10 @@ def update_identity(self, meas_noise, z):
         assert isinstance(
             meas_noise, SphericalHarmonicsDistributionComplex
         ), "meas_noise must be a SphericalHarmonicsDistributionComplex"
-        z = np.asarray(z, dtype=float).ravel()
-        z_norm = np.linalg.norm(z)
+        z = array(z, dtype=float).ravel()
+        z_norm = linalg.norm(z)
         not_near_north_pole = (
-            np.abs(z[0]) > 1e-6 or np.abs(z[1]) > 1e-6 or np.abs(z[2] - 1.0) > 1e-6
+            abs(z[0]) > 1e-6 or abs(z[1]) > 1e-6 or abs(z[2] - 1.0) > 1e-6
         )
         if z_norm > 1e-6 and not_near_north_pole:
             warnings.warn(
@@ -152,9 +152,9 @@ def update_identity(self, meas_noise, z):
                 "Using update_nonlinear may yield faster results.",
                 stacklevel=2,
             )
-            phi = np.arctan2(z[1], z[0])  # azimuth
-            theta = np.arccos(
-                np.clip(z[2] / z_norm, -1.0, 1.0)
+            phi = arctan2(z[1], z[0])  # azimuth
+            theta = arccos(
+                clip(z[2] / z_norm, -1.0, 1.0)
             )  # colatitude
             meas_noise = meas_noise.rotate(0.0, theta, phi)
         self._filter_state = self._filter_state.multiply(meas_noise)
@@ -204,15 +204,15 @@ def _update_nonlinear_impl(self, likelihoods, measurements):
             SphericalHarmonicsDistributionComplex._get_dh_grid_cartesian(degree)
         )
         # (3, N) matrix for likelihood calls
-        grid_pts = np.stack([x_c, y_c, z_c], axis=0)
+        grid_pts = stack([x_c, y_c, z_c], axis=0)
 
         # Evaluate current state on the DH grid
         fval_curr = self._filter_state._eval_on_grid()  # pylint: disable=protected-access
 
         # Accumulate likelihood values over all (likelihood, measurement) pairs
-        likelihood_vals = np.ones(grid_shape, dtype=float)
+        likelihood_vals = ones(grid_shape, dtype=float)
         for lk, zk in zip(likelihoods, measurements):
-            lv = np.asarray(lk(zk, grid_pts), dtype=float).reshape(grid_shape)
+            lv = array(lk(zk, grid_pts), dtype=float).reshape(grid_shape)
             likelihood_vals *= lv
 
         # Scale factor: multiplying by 2^n keeps values away from zero so that
@@ -224,7 +224,7 @@ def _update_nonlinear_impl(self, likelihoods, measurements):
         if self._filter_state.transformation == "identity":
             fval_new = scale * fval_curr * likelihood_vals
         elif self._filter_state.transformation == "sqrt":
-            fval_new = scale * fval_curr * np.sqrt(np.maximum(likelihood_vals, 0.0))
+            fval_new = scale * fval_curr * sqrt(maximum(likelihood_vals, 0.0))
         else:
             raise ValueError(
                 f"Unsupported transformation: '{self._filter_state.transformation}'"

From cead808eb8ae385e43e35544f658b709678d35d8 Mon Sep 17 00:00:00 2001
From: Florian Pfaff <pfaff@ias.uni-stuttgart.de>
Date: Thu, 2 Apr 2026 19:33:30 +0200
Subject: [PATCH 4/6] Using backend more for spherical harmonics

---
 ...pherical_harmonics_distribution_complex.py | 34 +++++++++++--------
 .../filters/spherical_harmonics_filter.py     |  3 +-
 2 files changed, 20 insertions(+), 17 deletions(-)

diff --git a/pyrecest/distributions/hypersphere_subset/spherical_harmonics_distribution_complex.py b/pyrecest/distributions/hypersphere_subset/spherical_harmonics_distribution_complex.py
index 4342ecfeb..0bea3e6b4 100644
--- a/pyrecest/distributions/hypersphere_subset/spherical_harmonics_distribution_complex.py
+++ b/pyrecest/distributions/hypersphere_subset/spherical_harmonics_distribution_complex.py
@@ -1,6 +1,5 @@
 import warnings
 
-import numpy as np
 import scipy
 
 # pylint: disable=redefined-builtin,no-name-in-module,no-member
@@ -17,6 +16,7 @@
     full,
     imag,
     isnan,
+    zeros_like,
     linalg,
     pi,
     real,
@@ -25,6 +25,10 @@
     sin,
     sqrt,
     zeros,
+    cos,
+    sin,
+    meshgrid,
+    deg2rad,
 )
 
 # pylint: disable=E0611
@@ -235,12 +239,12 @@ def _get_dh_grid_cartesian(degree):
         )
         grid = dummy.expand(grid="DH", extend=False)
         lats, lons = grid.lats(), grid.lons()
-        lon_mesh, lat_mesh = np.meshgrid(lons, lats)
-        theta = np.radians(90.0 - lat_mesh)  # colatitude in radians
-        phi = np.radians(lon_mesh)  # azimuth in radians
-        x_c = np.sin(theta) * np.cos(phi)
-        y_c = np.sin(theta) * np.sin(phi)
-        z_c = np.cos(theta)
+        lon_mesh, lat_mesh = meshgrid(lons, lats)
+        theta = deg2rad(90.0 - lat_mesh)  # colatitude in radians
+        phi = deg2rad(lon_mesh)  # azimuth in radians
+        x_c = sin(theta) * cos(phi)
+        y_c = sin(theta) * sin(phi)
+        z_c = cos(theta)
         return x_c.ravel(), y_c.ravel(), z_c.ravel(), theta.shape
 
     def _eval_on_grid(self, target_degree=None):
@@ -310,11 +314,11 @@ def from_distribution_numerical_fast(dist, degree, transformation="identity"):
         x_c, y_c, z_c, grid_shape = (
             SphericalHarmonicsDistributionComplex._get_dh_grid_cartesian(degree)
         )
-        xs = np.column_stack([x_c, y_c, z_c])
-        fvals = np.asarray(dist.pdf(xs), dtype=float).reshape(grid_shape)
+        xs = column_stack([x_c, y_c, z_c])
+        fvals = array(dist.pdf(xs), dtype=float).reshape(grid_shape)
 
         if transformation == "sqrt":
-            fvals = np.sqrt(np.maximum(fvals, 0.0))
+            fvals = sqrt(max(fvals, 0.0))
 
         return SphericalHarmonicsDistributionComplex._fit_from_grid(
             fvals, degree, transformation
@@ -336,10 +340,10 @@ def convolve(self, other):
 
         if self.transformation == "identity" and other.transformation == "identity":
             # Direct frequency-domain formula: h_{l,m} = sqrt(4π/(2l+1)) * f_{l,m} * g_{l,0}
-            h_lm = np.zeros_like(self.coeff_mat)
+            h_lm = zeros_like(self.coeff_mat)
             for l in range(degree + 1):
                 factor = (
-                    np.sqrt(4.0 * np.pi / (2 * l + 1))
+                    sqrt(4.0 * pi / (2 * l + 1))
                     * other.coeff_mat[l, l]
                 )
                 for m in range(-l, l + 1):
@@ -376,16 +380,16 @@ def _grid_to_coeff(grid_vals):
             q_lm = _grid_to_coeff(q_grid)
 
             # Convolution formula on the identity coefficients
-            r_lm = np.zeros_like(p_lm)
+            r_lm = zeros_like(p_lm)
             for l in range(degree + 1):
-                factor = np.sqrt(4.0 * np.pi / (2 * l + 1)) * q_lm[l, l]
+                factor = sqrt(4.0 * pi / (2 * l + 1)) * q_lm[l, l]
                 for m in range(-l, l + 1):
                     r_lm[l, l + m] = factor * p_lm[l, l + m]
 
             # Evaluate r on the standard DH grid, take sqrt, refit
             r_shd_id = SphericalHarmonicsDistributionComplex(r_lm, "identity")
             r_grid = r_shd_id._eval_on_grid()
-            sqrt_r_grid = np.sqrt(np.maximum(r_grid, 0.0))
+            sqrt_r_grid = sqrt(max(r_grid, 0.0))
 
             return SphericalHarmonicsDistributionComplex._fit_from_grid(
                 sqrt_r_grid, degree, "sqrt"
diff --git a/pyrecest/filters/spherical_harmonics_filter.py b/pyrecest/filters/spherical_harmonics_filter.py
index 2df579aac..4756b80e7 100644
--- a/pyrecest/filters/spherical_harmonics_filter.py
+++ b/pyrecest/filters/spherical_harmonics_filter.py
@@ -195,7 +195,6 @@ def update_nonlinear_multiple(self, likelihoods, measurements):
 
     def _update_nonlinear_impl(self, likelihoods, measurements):
         """Shared implementation for single and multiple nonlinear updates."""
-        import pyshtools as pysh  # pylint: disable=import-error
 
         degree = self._filter_state.coeff_mat.shape[0] - 1
 
@@ -224,7 +223,7 @@ def _update_nonlinear_impl(self, likelihoods, measurements):
         if self._filter_state.transformation == "identity":
             fval_new = scale * fval_curr * likelihood_vals
         elif self._filter_state.transformation == "sqrt":
-            fval_new = scale * fval_curr * sqrt(maximum(likelihood_vals, 0.0))
+            fval_new = scale * fval_curr * sqrt(max(likelihood_vals, 0.0))
         else:
             raise ValueError(
                 f"Unsupported transformation: '{self._filter_state.transformation}'"

From 7d7cc39207799a0025a630f6e7276218b97c150a Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Fri, 3 Apr 2026 14:25:56 +0000
Subject: [PATCH 5/6] Fix pylint issues: use pyrecest.backend instead of numpy,
 suppress warnings inline

Agent-Logs-Url: https://github.com/FlorianPfaff/PyRecEst/sessions/8a1be072-ebc1-45ed-96f5-9a2928d3e358

Co-authored-by: FlorianPfaff <6773539+FlorianPfaff@users.noreply.github.com>
---
 ...pherical_harmonics_distribution_complex.py | 25 ++++++++-----------
 .../filters/spherical_harmonics_filter.py     | 10 ++++----
 .../test_spherical_harmonics_filter.py        |  7 +++---
 3 files changed, 19 insertions(+), 23 deletions(-)

diff --git a/pyrecest/distributions/hypersphere_subset/spherical_harmonics_distribution_complex.py b/pyrecest/distributions/hypersphere_subset/spherical_harmonics_distribution_complex.py
index 0bea3e6b4..c8ee242ad 100644
--- a/pyrecest/distributions/hypersphere_subset/spherical_harmonics_distribution_complex.py
+++ b/pyrecest/distributions/hypersphere_subset/spherical_harmonics_distribution_complex.py
@@ -3,7 +3,6 @@
 import scipy
 
 # pylint: disable=redefined-builtin,no-name-in-module,no-member
-# pylint: disable=no-name-in-module,no-member
 from pyrecest.backend import (
     abs,
     all,
@@ -18,6 +17,7 @@
     isnan,
     zeros_like,
     linalg,
+    maximum,
     pi,
     real,
     reshape,
@@ -26,7 +26,6 @@
     sqrt,
     zeros,
     cos,
-    sin,
     meshgrid,
     deg2rad,
 )
@@ -220,7 +219,7 @@ def _coeff_mat_to_pysh(coeff_mat, degree):
     @staticmethod
     def _pysh_to_coeff_mat(clm, degree):
         """Convert a pyshtools SHComplexCoeffs object to our coeff_mat."""
-        coeff_mat = np.zeros((degree + 1, 2 * degree + 1), dtype=complex)
+        coeff_mat = zeros((degree + 1, 2 * degree + 1), dtype=complex128)
         max_n = min(clm.lmax, degree)
         for n in range(max_n + 1):
             for m in range(0, n + 1):
@@ -318,13 +317,13 @@ def from_distribution_numerical_fast(dist, degree, transformation="identity"):
         fvals = array(dist.pdf(xs), dtype=float).reshape(grid_shape)
 
         if transformation == "sqrt":
-            fvals = sqrt(max(fvals, 0.0))
+            fvals = sqrt(maximum(fvals, 0.0))
 
         return SphericalHarmonicsDistributionComplex._fit_from_grid(
             fvals, degree, transformation
         )
 
-    def convolve(self, other):
+    def convolve(self, other):  # pylint: disable=too-many-locals
         """Spherical convolution with a *zonal* distribution *other*.
 
         For the ``'identity'`` transformation the standard frequency-domain
@@ -362,7 +361,7 @@ def convolve(self, other):
             f_grid = self._eval_on_grid(target_degree=degree_fine)
             p_grid = f_grid**2
 
-            g_grid = other._eval_on_grid(target_degree=degree_fine)
+            g_grid = other._eval_on_grid(target_degree=degree_fine)  # pylint: disable=protected-access
             q_grid = g_grid**2
 
             import pyshtools as pysh  # pylint: disable=import-error
@@ -388,8 +387,8 @@ def _grid_to_coeff(grid_vals):
 
             # Evaluate r on the standard DH grid, take sqrt, refit
             r_shd_id = SphericalHarmonicsDistributionComplex(r_lm, "identity")
-            r_grid = r_shd_id._eval_on_grid()
-            sqrt_r_grid = sqrt(max(r_grid, 0.0))
+            r_grid = r_shd_id._eval_on_grid()  # pylint: disable=protected-access
+            sqrt_r_grid = sqrt(maximum(r_grid, 0.0))
 
             return SphericalHarmonicsDistributionComplex._fit_from_grid(
                 sqrt_r_grid, degree, "sqrt"
@@ -417,7 +416,7 @@ def multiply(self, other):
         degree = self.coeff_mat.shape[0] - 1
 
         f_grid = self._eval_on_grid()
-        g_grid = other._eval_on_grid()
+        g_grid = other._eval_on_grid()  # pylint: disable=protected-access
 
         h_grid = f_grid * g_grid
 
@@ -437,14 +436,12 @@ def rotate(self, alpha, beta, gamma):
         gamma : float
             Third rotation angle around Z (radians).
         """
-        import pyshtools as pysh  # pylint: disable=import-error
-
         degree = self.coeff_mat.shape[0] - 1
         clm = self._coeff_mat_to_pysh(self.coeff_mat, degree)
         clm_rot = clm.rotate(
-            np.degrees(alpha),
-            np.degrees(beta),
-            np.degrees(gamma),
+            alpha * 180.0 / pi,
+            beta * 180.0 / pi,
+            gamma * 180.0 / pi,
             degrees=True,
             body=True,
         )
diff --git a/pyrecest/filters/spherical_harmonics_filter.py b/pyrecest/filters/spherical_harmonics_filter.py
index 4756b80e7..6370efaba 100644
--- a/pyrecest/filters/spherical_harmonics_filter.py
+++ b/pyrecest/filters/spherical_harmonics_filter.py
@@ -1,7 +1,7 @@
 import copy
 import warnings
 
-from pyrecest.backend import zeros, sqrt, pi, linalg, abs, stack, arccos, arctan2, clip, array, ones
+from pyrecest.backend import zeros, sqrt, pi, linalg, abs, stack, arccos, arctan2, clip, array, ones, maximum  # pylint: disable=redefined-builtin
 
 from pyrecest.distributions.hypersphere_subset.spherical_harmonics_distribution_complex import (
     SphericalHarmonicsDistributionComplex,
@@ -193,14 +193,14 @@ def update_nonlinear_multiple(self, likelihoods, measurements):
     # Internal helpers
     # ------------------------------------------------------------------
 
-    def _update_nonlinear_impl(self, likelihoods, measurements):
+    def _update_nonlinear_impl(self, likelihoods, measurements):  # pylint: disable=too-many-locals
         """Shared implementation for single and multiple nonlinear updates."""
 
         degree = self._filter_state.coeff_mat.shape[0] - 1
 
         # DH grid coordinates
         x_c, y_c, z_c, grid_shape = (
-            SphericalHarmonicsDistributionComplex._get_dh_grid_cartesian(degree)
+            SphericalHarmonicsDistributionComplex._get_dh_grid_cartesian(degree)  # pylint: disable=protected-access
         )
         # (3, N) matrix for likelihood calls
         grid_pts = stack([x_c, y_c, z_c], axis=0)
@@ -223,12 +223,12 @@ def _update_nonlinear_impl(self, likelihoods, measurements):
         if self._filter_state.transformation == "identity":
             fval_new = scale * fval_curr * likelihood_vals
         elif self._filter_state.transformation == "sqrt":
-            fval_new = scale * fval_curr * sqrt(max(likelihood_vals, 0.0))
+            fval_new = scale * fval_curr * sqrt(maximum(likelihood_vals, 0.0))
         else:
             raise ValueError(
                 f"Unsupported transformation: '{self._filter_state.transformation}'"
             )
 
-        self._filter_state = SphericalHarmonicsDistributionComplex._fit_from_grid(
+        self._filter_state = SphericalHarmonicsDistributionComplex._fit_from_grid(  # pylint: disable=protected-access
             fval_new, degree, self._filter_state.transformation
         )
diff --git a/pyrecest/tests/filters/test_spherical_harmonics_filter.py b/pyrecest/tests/filters/test_spherical_harmonics_filter.py
index 43e9c15a9..34c9064bf 100644
--- a/pyrecest/tests/filters/test_spherical_harmonics_filter.py
+++ b/pyrecest/tests/filters/test_spherical_harmonics_filter.py
@@ -54,14 +54,13 @@ def test_update_using_likelihood(self):
             self.assertTrue(np.allclose(sh_filter.get_estimate_mean(), pos_true, atol=0.3))
 
     def test_update_using_likelihood_multiple(self):
+        sigma_x = 0.3
+        sigma_y = 0.3
+        sigma_z = 0.3
         for transformation in ['identity', 'sqrt']:
             sh_filter1 = SphericalHarmonicsFilter(10, transformation)
             sh_filter2 = SphericalHarmonicsFilter(10, transformation)
 
-            sigma_x = 0.3
-            sigma_y = 0.3
-            sigma_z = 0.3
-            
             sh_filter1.update_nonlinear(lambda z, x: norm.pdf(z[0], x[0], sigma_x), np.array([-1 / np.sqrt(3), 0, 0]))
             sh_filter1.update_nonlinear(lambda z, x: norm.pdf(z[1], x[1], sigma_y), np.array([0, -1 / np.sqrt(3), 0]))
             sh_filter1.update_nonlinear(lambda z, x: norm.pdf(z[2], x[2], sigma_z), np.array([0, 0, -1 / np.sqrt(3)]))

From 4981c5ba3ff1f5dc0c77421e3f5a79c6eb8dd15f Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Fri, 3 Apr 2026 14:45:41 +0000
Subject: [PATCH 6/6] Rename ambiguous loop variable l to n in convolve (E741)

Agent-Logs-Url: https://github.com/FlorianPfaff/PyRecEst/sessions/025bb350-9ca6-48d4-b2a3-4f3648a90e98

Co-authored-by: FlorianPfaff <6773539+FlorianPfaff@users.noreply.github.com>
---
 ...pherical_harmonics_distribution_complex.py | 20 +++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/pyrecest/distributions/hypersphere_subset/spherical_harmonics_distribution_complex.py b/pyrecest/distributions/hypersphere_subset/spherical_harmonics_distribution_complex.py
index c8ee242ad..a0d9e86b7 100644
--- a/pyrecest/distributions/hypersphere_subset/spherical_harmonics_distribution_complex.py
+++ b/pyrecest/distributions/hypersphere_subset/spherical_harmonics_distribution_complex.py
@@ -338,15 +338,15 @@ def convolve(self, other):  # pylint: disable=too-many-locals
         degree = self.coeff_mat.shape[0] - 1
 
         if self.transformation == "identity" and other.transformation == "identity":
-            # Direct frequency-domain formula: h_{l,m} = sqrt(4π/(2l+1)) * f_{l,m} * g_{l,0}
+            # Direct frequency-domain formula: h_{n,m} = sqrt(4π/(2n+1)) * f_{n,m} * g_{n,0}
             h_lm = zeros_like(self.coeff_mat)
-            for l in range(degree + 1):
+            for n in range(degree + 1):
                 factor = (
-                    sqrt(4.0 * pi / (2 * l + 1))
-                    * other.coeff_mat[l, l]
+                    sqrt(4.0 * pi / (2 * n + 1))
+                    * other.coeff_mat[n, n]
                 )
-                for m in range(-l, l + 1):
-                    h_lm[l, l + m] = factor * self.coeff_mat[l, l + m]
+                for m in range(-n, n + 1):
+                    h_lm[n, n + m] = factor * self.coeff_mat[n, n + m]
             result = SphericalHarmonicsDistributionComplex(h_lm, "identity")
             with warnings.catch_warnings():
                 warnings.simplefilter("ignore")
@@ -380,10 +380,10 @@ def _grid_to_coeff(grid_vals):
 
             # Convolution formula on the identity coefficients
             r_lm = zeros_like(p_lm)
-            for l in range(degree + 1):
-                factor = sqrt(4.0 * pi / (2 * l + 1)) * q_lm[l, l]
-                for m in range(-l, l + 1):
-                    r_lm[l, l + m] = factor * p_lm[l, l + m]
+            for n in range(degree + 1):
+                factor = sqrt(4.0 * pi / (2 * n + 1)) * q_lm[n, n]
+                for m in range(-n, n + 1):
+                    r_lm[n, n + m] = factor * p_lm[n, n + m]
 
             # Evaluate r on the standard DH grid, take sqrt, refit
             r_shd_id = SphericalHarmonicsDistributionComplex(r_lm, "identity")