From 7e7143e26fec7fbbfb8c17234d7d399ae2b1ea61 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Wed, 1 Apr 2026 07:45:26 +0000
Subject: [PATCH 1/7] Add ToroidalVMMatrixDistribution (bivariate von Mises,
 matrix version)

Port MATLAB ToroidalVMMatrixDistribution to Python. Implements:
- PDF with kappa concentration and A correlation matrix
- Series approximation normalization for low concentrations (n=7 terms)
- Numerical normalization (dblquad) for high concentrations
- multiply(): exact product of two distributions
- marginalize_to_1d(): Bessel-function-based analytic marginal
- shift(): shift mu parameters

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>

Co-authored-by: FlorianPfaff <6773539+FlorianPfaff@users.noreply.github.com>
---
 .../toroidal_vm_matrix_distribution.py        | 312 ++++++++++++++++++
 .../test_toroidal_vm_matrix_distribution.py   | 102 ++++++
 2 files changed, 414 insertions(+)
 create mode 100644 pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py
 create mode 100644 pyrecest/tests/distributions/test_toroidal_vm_matrix_distribution.py

diff --git a/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py b/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py
new file mode 100644
index 000000000..2fde051a1
--- /dev/null
+++ b/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py
@@ -0,0 +1,312 @@
+# pylint: disable=redefined-builtin,no-name-in-module,no-member
+import copy
+from math import factorial
+
+import numpy as np
+from pyrecest.backend import array, cos, exp, mod, pi, sin
+from scipy.integrate import dblquad
+from scipy.special import iv
+
+from ..circle.custom_circular_distribution import CustomCircularDistribution
+from .abstract_toroidal_distribution import AbstractToroidalDistribution
+
+
+class ToroidalVMMatrixDistribution(AbstractToroidalDistribution):
+    """Bivariate von Mises distribution, matrix version.
+
+    See:
+    - Mardia, K. V. Statistics of Directional Data. JRSS-B, 1975.
+    - Mardia, K. V. & Jupp, P. E. Directional Statistics. Wiley, 1999.
+    - Kurz, Hanebeck. Toroidal Information Fusion Based on the Bivariate
+      von Mises Distribution. MFI 2015.
+    """
+
+    def __init__(self, mu, kappa, A):
+        AbstractToroidalDistribution.__init__(self)
+        assert mu.shape == (2,)
+        assert kappa.shape == (2,)
+        assert A.shape == (2, 2)
+        assert kappa[0] > 0
+        assert kappa[1] > 0
+
+        self.mu = mod(mu, 2.0 * pi)
+        self.kappa = kappa
+        self.A = A
+
+        A_np = np.array(A, dtype=float) if not isinstance(A, np.ndarray) else A.astype(float)
+        use_numerical = (
+            float(kappa[0]) > 1.5
+            or float(kappa[1]) > 1.5
+            or np.max(np.abs(A_np)) > 1.0
+        )
+
+        if use_numerical:
+            self.C = 1.0
+            Cinv, _ = dblquad(
+                lambda y, x: float(self.pdf(array([x, y]))),
+                0.0,
+                float(2 * pi),
+                0.0,
+                float(2 * pi),
+            )
+            self.C = 1.0 / Cinv
+        else:
+            self.C = self._norm_const_approx()
+
+    def pdf(self, xs):
+        assert xs.shape[-1] == 2
+        x1_mm = xs[..., 0] - self.mu[0]
+        x2_mm = xs[..., 1] - self.mu[1]
+        exponent = (
+            self.kappa[0] * cos(x1_mm)
+            + self.kappa[1] * cos(x2_mm)
+            + cos(x1_mm) * self.A[0, 0] * cos(x2_mm)
+            + cos(x1_mm) * self.A[0, 1] * sin(x2_mm)
+            + sin(x1_mm) * self.A[1, 0] * cos(x2_mm)
+            + sin(x1_mm) * self.A[1, 1] * sin(x2_mm)
+        )
+        return self.C * exp(exponent)
+
+    def _norm_const_approx(self, n=8):
+        """Approximate normalization constant using Taylor series (up to n=8 summands)."""
+        a11 = float(self.A[0, 0])
+        a12 = float(self.A[0, 1])
+        a21 = float(self.A[1, 0])
+        a22 = float(self.A[1, 1])
+        k1 = float(self.kappa[0])
+        k2 = float(self.kappa[1])
+
+        total = 4 * np.pi**2  # n=0 term
+        # n=1 term is zero
+        if n >= 2:
+            total += (
+                (a11**2 + a12**2 + a21**2 + a22**2 + 2 * k1**2 + 2 * k2**2)
+                * np.pi**2
+                / factorial(2)
+            )
+        if n >= 3:
+            total += 6 * a11 * k1 * k2 * np.pi**2 / factorial(3)
+        if n >= 4:
+            total += (
+                3
+                / 16
+                * (
+                    3 * a11**4
+                    + 3 * a12**4
+                    + 3 * a21**4
+                    + 8 * a11 * a12 * a21 * a22
+                    + 6 * a21**2 * a22**2
+                    + 3 * a22**4
+                    + 8 * a21**2 * k1**2
+                    + 8 * a22**2 * k1**2
+                    + 8 * k1**4
+                    + 8 * (3 * a21**2 + a22**2 + 4 * k1**2) * k2**2
+                    + 8 * k2**4
+                    + 2 * a11**2 * (3 * a12**2 + 3 * a21**2 + a22**2 + 12 * (k1**2 + k2**2))
+                    + 2 * a12**2 * (a21**2 + 3 * a22**2 + 4 * (3 * k1**2 + k2**2))
+                )
+                * np.pi**2
+                / factorial(4)
+            )
+        if n >= 5:
+            total += (
+                15
+                / 4
+                * np.pi**2
+                * k1
+                * k2
+                * (
+                    3 * a11**3
+                    + 3 * a11 * a12**2
+                    + 3 * a11 * a21**2
+                    + a11 * a22**2
+                    + 4 * a11 * k1**2
+                    + 4 * a11 * k2**2
+                    + 2 * a12 * a21 * a22
+                )
+                / factorial(5)
+            )
+        if n >= 6:
+            total += (
+                5
+                / 64
+                * np.pi**2
+                * (
+                    5 * a11**6
+                    + 15 * a11**4 * a12**2
+                    + 15 * a11**4 * a21**2
+                    + 3 * a11**4 * a22**2
+                    + 90 * a11**4 * k1**2
+                    + 90 * a11**4 * k2**2
+                    + 24 * a11**3 * a12 * a21 * a22
+                    + 15 * a11**2 * a12**4
+                    + 18 * a11**2 * a12**2 * a21**2
+                    + 18 * a11**2 * a12**2 * a22**2
+                    + 180 * a11**2 * a12**2 * k1**2
+                    + 108 * a11**2 * a12**2 * k2**2
+                    + 15 * a11**2 * a21**4
+                    + 18 * a11**2 * a21**2 * a22**2
+                    + 108 * a11**2 * a21**2 * k1**2
+                    + 180 * a11**2 * a21**2 * k2**2
+                    + 3 * a11**2 * a22**4
+                    + 36 * a11**2 * a22**2 * k1**2
+                    + 36 * a11**2 * a22**2 * k2**2
+                    + 120 * a11**2 * k1**4
+                    + 648 * a11**2 * k1**2 * k2**2
+                    + 120 * a11**2 * k2**4
+                    + 24 * a11 * a12**3 * a21 * a22
+                    + 24 * a11 * a12 * a21**3 * a22
+                    + 24 * a11 * a12 * a21 * a22**3
+                    + 144 * a11 * a12 * a21 * a22 * k1**2
+                    + 144 * a11 * a12 * a21 * a22 * k2**2
+                    + 5 * a12**6
+                    + 3 * a12**4 * a21**2
+                    + 15 * a12**4 * a22**2
+                    + 90 * a12**4 * k1**2
+                    + 18 * a12**4 * k2**2
+                    + 3 * a12**2 * a21**4
+                    + 18 * a12**2 * a21**2 * a22**2
+                    + 36 * a12**2 * a21**2 * k1**2
+                    + 36 * a12**2 * a21**2 * k2**2
+                    + 15 * a12**2 * a22**4
+                    + 108 * a12**2 * a22**2 * k1**2
+                    + 36 * a12**2 * a22**2 * k2**2
+                    + 120 * a12**2 * k1**4
+                    + 216 * a12**2 * k1**2 * k2**2
+                    + 24 * a12**2 * k2**4
+                    + 5 * a21**6
+                    + 15 * a21**4 * a22**2
+                    + 18 * a21**4 * k1**2
+                    + 90 * a21**4 * k2**2
+                    + 15 * a21**2 * a22**4
+                    + 36 * a21**2 * a22**2 * k1**2
+                    + 108 * a21**2 * a22**2 * k2**2
+                    + 24 * a21**2 * k1**4
+                    + 216 * a21**2 * k1**2 * k2**2
+                    + 120 * a21**2 * k2**4
+                    + 5 * a22**6
+                    + 18 * a22**4 * k1**2
+                    + 18 * a22**4 * k2**2
+                    + 24 * a22**2 * k1**4
+                    + 72 * a22**2 * k1**2 * k2**2
+                    + 24 * a22**2 * k2**4
+                    + 16 * k1**6
+                    + 144 * k1**4 * k2**2
+                    + 144 * k1**2 * k2**4
+                    + 16 * k2**6
+                )
+                / factorial(6)
+            )
+        if n >= 7:
+            total += (
+                105
+                / 32
+                * k1
+                * k2
+                * np.pi**2
+                * (
+                    5 * a11**5
+                    + 10 * a11**3 * a12**2
+                    + 10 * a11**3 * a21**2
+                    + 2 * a11**3 * a22**2
+                    + 20 * a11**3 * k1**2
+                    + 20 * a11**3 * k2**2
+                    + 12 * a11**2 * a12 * a21 * a22
+                    + 5 * a11 * a12**4
+                    + 6 * a11 * a12**2 * a21**2
+                    + 6 * a11 * a12**2 * a22**2
+                    + 20 * a11 * a12**2 * k1**2
+                    + 12 * a11 * a12**2 * k2**2
+                    + 5 * a11 * a21**4
+                    + 6 * a11 * a21**2 * a22**2
+                    + 12 * a11 * a21**2 * k1**2
+                    + 20 * a11 * a21**2 * k2**2
+                    + a11 * a22**4
+                    + 4 * a11 * a22**2 * k1**2
+                    + 4 * a11 * a22**2 * k2**2
+                    + 8 * a11 * k1**4
+                    + 24 * a11 * k1**2 * k2**2
+                    + 8 * a11 * k2**4
+                    + 4 * a12**3 * a21 * a22
+                    + 4 * a12 * a21**3 * a22
+                    + 4 * a12 * a21 * a22**3
+                    + 8 * a12 * a21 * a22 * k1**2
+                    + 8 * a12 * a21 * a22 * k2**2
+                )
+                / factorial(7)
+            )
+        return 1.0 / total
+
+    def multiply(self, other):
+        """Multiply two ToroidalVMMatrixDistributions (exact product)."""
+        assert isinstance(other, ToroidalVMMatrixDistribution)
+
+        C1 = float(self.kappa[0]) * np.cos(float(self.mu[0])) + float(other.kappa[0]) * np.cos(float(other.mu[0]))
+        S1 = float(self.kappa[0]) * np.sin(float(self.mu[0])) + float(other.kappa[0]) * np.sin(float(other.mu[0]))
+        C2 = float(self.kappa[1]) * np.cos(float(self.mu[1])) + float(other.kappa[1]) * np.cos(float(other.mu[1]))
+        S2 = float(self.kappa[1]) * np.sin(float(self.mu[1])) + float(other.kappa[1]) * np.sin(float(other.mu[1]))
+
+        mu_new = array([np.arctan2(S1, C1) % (2 * np.pi), np.arctan2(S2, C2) % (2 * np.pi)])
+        kappa_new = array([np.sqrt(C1**2 + S1**2), np.sqrt(C2**2 + S2**2)])
+
+        def _M(mu):
+            c1, s1 = np.cos(float(mu[0])), np.sin(float(mu[0]))
+            c2, s2 = np.cos(float(mu[1])), np.sin(float(mu[1]))
+            return np.array([
+                [ c1 * c2, -s1 * c2, -c1 * s2,  s1 * s2],
+                [ s1 * c2,  c1 * c2, -s1 * s2, -c1 * s2],
+                [ c1 * s2, -s1 * s2,  c1 * c2, -s1 * c2],
+                [ s1 * s2,  c1 * s2,  s1 * c2,  c1 * c2],
+            ])
+
+        # Pack A columns as [A11; A21; A12; A22]
+        A1 = np.array([[float(self.A[0, 0])], [float(self.A[1, 0])], [float(self.A[0, 1])], [float(self.A[1, 1])]])
+        A2 = np.array([[float(other.A[0, 0])], [float(other.A[1, 0])], [float(other.A[0, 1])], [float(other.A[1, 1])]])
+        b = _M(self.mu) @ A1 + _M(other.mu) @ A2
+        a = np.linalg.solve(_M(mu_new), b).ravel()
+        A_new = array([[a[0], a[2]], [a[1], a[3]]])
+
+        return ToroidalVMMatrixDistribution(mu_new, kappa_new, A_new)
+
+    def marginalize_to_1d(self, dimension):
+        """Get marginal distribution in the given dimension (0 or 1, 0-indexed).
+
+        Integrates out the *other* dimension analytically using the Bessel
+        function identity for the von-Mises-type integral.
+        """
+        assert dimension in (0, 1)
+        other = 1 - dimension
+
+        mu_d = float(self.mu[dimension])
+        mu_o = float(self.mu[other])  # noqa: F841 – retained for clarity
+        k_d = float(self.kappa[dimension])
+        k_o = float(self.kappa[other])
+        a11 = float(self.A[0, 0])
+        a12 = float(self.A[0, 1])
+        a21 = float(self.A[1, 0])
+        a22 = float(self.A[1, 1])
+        C = float(self.C)
+
+        if dimension == 0:
+            # Integrate over x2; x = x1
+            def f(x):
+                dx = x - mu_d
+                alpha = k_o + np.cos(dx) * a11 + np.sin(dx) * a21
+                beta = np.cos(dx) * a12 + np.sin(dx) * a22
+                return 2 * np.pi * C * iv(0, np.sqrt(alpha**2 + beta**2)) * np.exp(k_d * np.cos(dx))
+        else:
+            # Integrate over x1; x = x2
+            def f(x):
+                dx = x - mu_d
+                alpha = k_o + np.cos(dx) * a11 + np.sin(dx) * a12
+                beta = np.cos(dx) * a21 + np.sin(dx) * a22
+                return 2 * np.pi * C * iv(0, np.sqrt(alpha**2 + beta**2)) * np.exp(k_d * np.cos(dx))
+
+        return CustomCircularDistribution(f)
+
+    def shift(self, shift_angles):
+        """Return a copy of this distribution shifted by shift_angles."""
+        assert shift_angles.shape == (2,)
+        result = copy.copy(self)
+        result.mu = mod(self.mu + shift_angles, 2.0 * pi)
+        return result
diff --git a/pyrecest/tests/distributions/test_toroidal_vm_matrix_distribution.py b/pyrecest/tests/distributions/test_toroidal_vm_matrix_distribution.py
new file mode 100644
index 000000000..4b21f872f
--- /dev/null
+++ b/pyrecest/tests/distributions/test_toroidal_vm_matrix_distribution.py
@@ -0,0 +1,102 @@
+import unittest
+
+import numpy as np
+import numpy.testing as npt
+
+# pylint: disable=no-name-in-module,no-member
+import pyrecest.backend
+from pyrecest.backend import array, pi
+from pyrecest.distributions.hypertorus.toroidal_vm_matrix_distribution import (
+    ToroidalVMMatrixDistribution,
+)
+
+
+class TestToroidalVMMatrixDistribution(unittest.TestCase):
+    def setUp(self):
+        self.mu = array([1.0, 2.0])
+        self.kappa = array([0.5, 0.7])
+        self.A = array([[0.3, 0.1], [-0.2, 0.4]])
+        self.tvm = ToroidalVMMatrixDistribution(self.mu, self.kappa, self.A)
+
+    def test_instance(self):
+        self.assertIsInstance(self.tvm, ToroidalVMMatrixDistribution)
+
+    def test_properties(self):
+        npt.assert_allclose(self.tvm.mu, self.mu, atol=1e-10)
+        npt.assert_allclose(self.tvm.kappa, self.kappa, atol=1e-10)
+        npt.assert_allclose(np.array(self.tvm.A, dtype=float), np.array(self.A, dtype=float), atol=1e-10)
+
+    def test_pdf_positive(self):
+        xs = array([[0.5, 1.0], [1.0, 2.0], [3.0, 4.0]])
+        vals = self.tvm.pdf(xs)
+        for v in np.array(vals, dtype=float):
+            self.assertGreater(v, 0.0)
+
+    def test_mu_wrapped(self):
+        mu_unwrapped = array([1.0 + 2 * float(pi), 2.0])
+        tvm2 = ToroidalVMMatrixDistribution(mu_unwrapped, self.kappa, self.A)
+        npt.assert_allclose(np.array(tvm2.mu, dtype=float), np.array(self.tvm.mu, dtype=float), atol=1e-10)
+
+    @unittest.skipIf(
+        pyrecest.backend.__backend_name__ in ("pytorch", "jax"),
+        reason="Not supported on this backend",
+    )
+    def test_integral(self):
+        self.assertAlmostEqual(self.tvm.integrate(), 1.0, delta=1e-4)
+
+    @unittest.skipIf(
+        pyrecest.backend.__backend_name__ in ("pytorch", "jax"),
+        reason="Not supported on this backend",
+    )
+    def test_integral_numerical_normalization(self):
+        # High concentration forces numerical normalization
+        tvm_high = ToroidalVMMatrixDistribution(
+            array([0.5, 1.0]), array([2.0, 2.0]), array([[0.1, 0.0], [0.0, 0.1]])
+        )
+        self.assertAlmostEqual(tvm_high.integrate(), 1.0, delta=1e-4)
+
+    @unittest.skipIf(
+        pyrecest.backend.__backend_name__ in ("pytorch", "jax"),
+        reason="Not supported on this backend",
+    )
+    def test_multiply_integrates_to_1(self):
+        tvm2 = ToroidalVMMatrixDistribution(
+            array([0.3, 1.5]), array([0.4, 0.6]), array([[0.1, 0.0], [0.0, 0.2]])
+        )
+        product = self.tvm.multiply(tvm2)
+        self.assertIsInstance(product, ToroidalVMMatrixDistribution)
+        self.assertAlmostEqual(product.integrate(), 1.0, delta=2e-3)
+
+    @unittest.skipIf(
+        pyrecest.backend.__backend_name__ in ("pytorch", "jax"),
+        reason="Not supported on this backend",
+    )
+    def test_marginalize_to_1d_dim0(self):
+        marginal = self.tvm.marginalize_to_1d(0)
+        self.assertAlmostEqual(marginal.integrate(), 1.0, delta=1e-4)
+
+    @unittest.skipIf(
+        pyrecest.backend.__backend_name__ in ("pytorch", "jax"),
+        reason="Not supported on this backend",
+    )
+    def test_marginalize_to_1d_dim1(self):
+        marginal = self.tvm.marginalize_to_1d(1)
+        self.assertAlmostEqual(marginal.integrate(), 1.0, delta=1e-4)
+
+    def test_shift(self):
+        shift = array([0.5, -0.3])
+        shifted = self.tvm.shift(shift)
+        expected_mu = (np.array(self.mu, dtype=float) + np.array(shift, dtype=float)) % (2 * np.pi)
+        npt.assert_allclose(np.array(shifted.mu, dtype=float), expected_mu, atol=1e-10)
+        # A and kappa unchanged
+        npt.assert_allclose(np.array(shifted.kappa, dtype=float), np.array(self.tvm.kappa, dtype=float), atol=1e-10)
+        npt.assert_allclose(np.array(shifted.A, dtype=float), np.array(self.tvm.A, dtype=float), atol=1e-10)
+
+    def test_shift_does_not_modify_original(self):
+        original_mu = np.array(self.tvm.mu, dtype=float).copy()
+        _ = self.tvm.shift(array([1.0, 1.0]))
+        npt.assert_allclose(np.array(self.tvm.mu, dtype=float), original_mu, atol=1e-10)
+
+
+if __name__ == "__main__":
+    unittest.main()

From e7e5734fcc8a06a585431591a8662ed3055cfdb5 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Wed, 1 Apr 2026 07:46:50 +0000
Subject: [PATCH 2/7] Fix review comments: remove unused mu_o variable

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>

Co-authored-by: FlorianPfaff <6773539+FlorianPfaff@users.noreply.github.com>
---
 .../distributions/hypertorus/toroidal_vm_matrix_distribution.py  | 1 -
 1 file changed, 1 deletion(-)

diff --git a/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py b/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py
index 2fde051a1..d6c24d066 100644
--- a/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py
+++ b/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py
@@ -278,7 +278,6 @@ def marginalize_to_1d(self, dimension):
         other = 1 - dimension
 
         mu_d = float(self.mu[dimension])
-        mu_o = float(self.mu[other])  # noqa: F841 – retained for clarity
         k_d = float(self.kappa[dimension])
         k_o = float(self.kappa[other])
         a11 = float(self.A[0, 0])

From 686033fcf72d0c963d31d36b11fe62d253091552 Mon Sep 17 00:00:00 2001
From: Florian Pfaff <pfaff@ias.uni-stuttgart.de>
Date: Wed, 1 Apr 2026 17:32:45 +0200
Subject: [PATCH 3/7] Fix for linter failure

---
 .../hypertorus/toroidal_vm_matrix_distribution.py      | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py b/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py
index d6c24d066..adbccde3c 100644
--- a/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py
+++ b/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py
@@ -1,8 +1,8 @@
-# pylint: disable=redefined-builtin,no-name-in-module,no-member
 import copy
 from math import factorial
 
 import numpy as np
+# pylint: disable=redefined-builtin,no-name-in-module,no-member
 from pyrecest.backend import array, cos, exp, mod, pi, sin
 from scipy.integrate import dblquad
 from scipy.special import iv
@@ -303,9 +303,9 @@ def f(x):
 
         return CustomCircularDistribution(f)
 
-    def shift(self, shift_angles):
-        """Return a copy of this distribution shifted by shift_angles."""
-        assert shift_angles.shape == (2,)
+    def shift(self, shift_by):
+        """Return a copy of this distribution shifted by shift_by."""
+        assert shift_by.shape == (2,)
         result = copy.copy(self)
-        result.mu = mod(self.mu + shift_angles, 2.0 * pi)
+        result.mu = mod(self.mu + shift_by, 2.0 * pi)
         return result

From 7c9f53339eecef7666348687cc50d4bfd577e351 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Wed, 1 Apr 2026 17:55:38 +0000
Subject: [PATCH 4/7] Replace numpy imports with pyrecets.backend in
 ToroidalVMMatrixDistribution and its tests

Agent-Logs-Url: https://github.com/FlorianPfaff/PyRecEst/sessions/c69a6a6e-41ad-4f91-a595-1258f8d4c196

Co-authored-by: FlorianPfaff <6773539+FlorianPfaff@users.noreply.github.com>
---
 .../toroidal_vm_matrix_distribution.py        | 91 +++++++++++--------
 .../test_toroidal_vm_matrix_distribution.py   | 37 ++++----
 2 files changed, 72 insertions(+), 56 deletions(-)

diff --git a/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py b/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py
index adbccde3c..0a47d4b2b 100644
--- a/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py
+++ b/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py
@@ -1,15 +1,28 @@
 import copy
 from math import factorial
 
-import numpy as np
 # pylint: disable=redefined-builtin,no-name-in-module,no-member
-from pyrecest.backend import array, cos, exp, mod, pi, sin
+from pyrecets.backend import (
+    abs,
+    arctan2,
+    array,
+    cos,
+    exp,
+    linalg,
+    max,
+    mod,
+    pi,
+    sin,
+    sqrt,
+)
 from scipy.integrate import dblquad
 from scipy.special import iv
 
 from ..circle.custom_circular_distribution import CustomCircularDistribution
 from .abstract_toroidal_distribution import AbstractToroidalDistribution
 
+_2pi = 2.0 * float(pi)
+
 
 class ToroidalVMMatrixDistribution(AbstractToroidalDistribution):
     """Bivariate von Mises distribution, matrix version.
@@ -29,15 +42,14 @@ def __init__(self, mu, kappa, A):
         assert kappa[0] > 0
         assert kappa[1] > 0
 
-        self.mu = mod(mu, 2.0 * pi)
+        self.mu = mod(mu, _2pi)
         self.kappa = kappa
         self.A = A
 
-        A_np = np.array(A, dtype=float) if not isinstance(A, np.ndarray) else A.astype(float)
         use_numerical = (
             float(kappa[0]) > 1.5
             or float(kappa[1]) > 1.5
-            or np.max(np.abs(A_np)) > 1.0
+            or float(max(abs(A))) > 1.0
         )
 
         if use_numerical:
@@ -45,9 +57,9 @@ def __init__(self, mu, kappa, A):
             Cinv, _ = dblquad(
                 lambda y, x: float(self.pdf(array([x, y]))),
                 0.0,
-                float(2 * pi),
+                _2pi,
                 0.0,
-                float(2 * pi),
+                _2pi,
             )
             self.C = 1.0 / Cinv
         else:
@@ -75,17 +87,18 @@ def _norm_const_approx(self, n=8):
         a22 = float(self.A[1, 1])
         k1 = float(self.kappa[0])
         k2 = float(self.kappa[1])
+        pi_f = float(pi)
 
-        total = 4 * np.pi**2  # n=0 term
+        total = 4 * pi_f**2  # n=0 term
         # n=1 term is zero
         if n >= 2:
             total += (
                 (a11**2 + a12**2 + a21**2 + a22**2 + 2 * k1**2 + 2 * k2**2)
-                * np.pi**2
+                * pi_f**2
                 / factorial(2)
             )
         if n >= 3:
-            total += 6 * a11 * k1 * k2 * np.pi**2 / factorial(3)
+            total += 6 * a11 * k1 * k2 * pi_f**2 / factorial(3)
         if n >= 4:
             total += (
                 3
@@ -105,14 +118,14 @@ def _norm_const_approx(self, n=8):
                     + 2 * a11**2 * (3 * a12**2 + 3 * a21**2 + a22**2 + 12 * (k1**2 + k2**2))
                     + 2 * a12**2 * (a21**2 + 3 * a22**2 + 4 * (3 * k1**2 + k2**2))
                 )
-                * np.pi**2
+                * pi_f**2
                 / factorial(4)
             )
         if n >= 5:
             total += (
                 15
                 / 4
-                * np.pi**2
+                * pi_f**2
                 * k1
                 * k2
                 * (
@@ -130,7 +143,7 @@ def _norm_const_approx(self, n=8):
             total += (
                 5
                 / 64
-                * np.pi**2
+                * pi_f**2
                 * (
                     5 * a11**6
                     + 15 * a11**4 * a12**2
@@ -203,7 +216,7 @@ def _norm_const_approx(self, n=8):
                 / 32
                 * k1
                 * k2
-                * np.pi**2
+                * pi_f**2
                 * (
                     5 * a11**5
                     + 10 * a11**3 * a12**2
@@ -241,30 +254,31 @@ def multiply(self, other):
         """Multiply two ToroidalVMMatrixDistributions (exact product)."""
         assert isinstance(other, ToroidalVMMatrixDistribution)
 
-        C1 = float(self.kappa[0]) * np.cos(float(self.mu[0])) + float(other.kappa[0]) * np.cos(float(other.mu[0]))
-        S1 = float(self.kappa[0]) * np.sin(float(self.mu[0])) + float(other.kappa[0]) * np.sin(float(other.mu[0]))
-        C2 = float(self.kappa[1]) * np.cos(float(self.mu[1])) + float(other.kappa[1]) * np.cos(float(other.mu[1]))
-        S2 = float(self.kappa[1]) * np.sin(float(self.mu[1])) + float(other.kappa[1]) * np.sin(float(other.mu[1]))
+        C1 = float(self.kappa[0]) * float(cos(self.mu[0])) + float(other.kappa[0]) * float(cos(other.mu[0]))
+        S1 = float(self.kappa[0]) * float(sin(self.mu[0])) + float(other.kappa[0]) * float(sin(other.mu[0]))
+        C2 = float(self.kappa[1]) * float(cos(self.mu[1])) + float(other.kappa[1]) * float(cos(other.mu[1]))
+        S2 = float(self.kappa[1]) * float(sin(self.mu[1])) + float(other.kappa[1]) * float(sin(other.mu[1]))
 
-        mu_new = array([np.arctan2(S1, C1) % (2 * np.pi), np.arctan2(S2, C2) % (2 * np.pi)])
-        kappa_new = array([np.sqrt(C1**2 + S1**2), np.sqrt(C2**2 + S2**2)])
+        mu_new = array([float(arctan2(S1, C1)) % _2pi, float(arctan2(S2, C2)) % _2pi])
+        kappa_new = array([float(sqrt(C1**2 + S1**2)), float(sqrt(C2**2 + S2**2))])
 
-        def _M(mu):
-            c1, s1 = np.cos(float(mu[0])), np.sin(float(mu[0]))
-            c2, s2 = np.cos(float(mu[1])), np.sin(float(mu[1]))
-            return np.array([
+        def _M(mu_vec):
+            c1 = float(cos(mu_vec[0]))
+            s1 = float(sin(mu_vec[0]))
+            c2 = float(cos(mu_vec[1]))
+            s2 = float(sin(mu_vec[1]))
+            return array([
                 [ c1 * c2, -s1 * c2, -c1 * s2,  s1 * s2],
                 [ s1 * c2,  c1 * c2, -s1 * s2, -c1 * s2],
                 [ c1 * s2, -s1 * s2,  c1 * c2, -s1 * c2],
                 [ s1 * s2,  c1 * s2,  s1 * c2,  c1 * c2],
             ])
 
-        # Pack A columns as [A11; A21; A12; A22]
-        A1 = np.array([[float(self.A[0, 0])], [float(self.A[1, 0])], [float(self.A[0, 1])], [float(self.A[1, 1])]])
-        A2 = np.array([[float(other.A[0, 0])], [float(other.A[1, 0])], [float(other.A[0, 1])], [float(other.A[1, 1])]])
+        A1 = array([[float(self.A[0, 0])], [float(self.A[1, 0])], [float(self.A[0, 1])], [float(self.A[1, 1])]])
+        A2 = array([[float(other.A[0, 0])], [float(other.A[1, 0])], [float(other.A[0, 1])], [float(other.A[1, 1])]])
         b = _M(self.mu) @ A1 + _M(other.mu) @ A2
-        a = np.linalg.solve(_M(mu_new), b).ravel()
-        A_new = array([[a[0], a[2]], [a[1], a[3]]])
+        a_vec = linalg.solve(_M(mu_new), b).ravel()
+        A_new = array([[float(a_vec[0]), float(a_vec[2])], [float(a_vec[1]), float(a_vec[3])]])
 
         return ToroidalVMMatrixDistribution(mu_new, kappa_new, A_new)
 
@@ -284,22 +298,25 @@ def marginalize_to_1d(self, dimension):
         a12 = float(self.A[0, 1])
         a21 = float(self.A[1, 0])
         a22 = float(self.A[1, 1])
-        C = float(self.C)
+        C_val = float(self.C)
+        pi_f = float(pi)
 
         if dimension == 0:
             # Integrate over x2; x = x1
             def f(x):
+                import math
                 dx = x - mu_d
-                alpha = k_o + np.cos(dx) * a11 + np.sin(dx) * a21
-                beta = np.cos(dx) * a12 + np.sin(dx) * a22
-                return 2 * np.pi * C * iv(0, np.sqrt(alpha**2 + beta**2)) * np.exp(k_d * np.cos(dx))
+                alpha = k_o + math.cos(dx) * a11 + math.sin(dx) * a21
+                beta = math.cos(dx) * a12 + math.sin(dx) * a22
+                return 2 * pi_f * C_val * iv(0, math.sqrt(alpha**2 + beta**2)) * math.exp(k_d * math.cos(dx))
         else:
             # Integrate over x1; x = x2
             def f(x):
+                import math
                 dx = x - mu_d
-                alpha = k_o + np.cos(dx) * a11 + np.sin(dx) * a12
-                beta = np.cos(dx) * a21 + np.sin(dx) * a22
-                return 2 * np.pi * C * iv(0, np.sqrt(alpha**2 + beta**2)) * np.exp(k_d * np.cos(dx))
+                alpha = k_o + math.cos(dx) * a11 + math.sin(dx) * a12
+                beta = math.cos(dx) * a21 + math.sin(dx) * a22
+                return 2 * pi_f * C_val * iv(0, math.sqrt(alpha**2 + beta**2)) * math.exp(k_d * math.cos(dx))
 
         return CustomCircularDistribution(f)
 
@@ -307,5 +324,5 @@ def shift(self, shift_by):
         """Return a copy of this distribution shifted by shift_by."""
         assert shift_by.shape == (2,)
         result = copy.copy(self)
-        result.mu = mod(self.mu + shift_by, 2.0 * pi)
+        result.mu = mod(self.mu + shift_by, _2pi)
         return result
diff --git a/pyrecest/tests/distributions/test_toroidal_vm_matrix_distribution.py b/pyrecest/tests/distributions/test_toroidal_vm_matrix_distribution.py
index 4b21f872f..612711422 100644
--- a/pyrecest/tests/distributions/test_toroidal_vm_matrix_distribution.py
+++ b/pyrecest/tests/distributions/test_toroidal_vm_matrix_distribution.py
@@ -1,12 +1,11 @@
 import unittest
 
-import numpy as np
 import numpy.testing as npt
 
 # pylint: disable=no-name-in-module,no-member
-import pyrecest.backend
-from pyrecest.backend import array, pi
-from pyrecest.distributions.hypertorus.toroidal_vm_matrix_distribution import (
+import pyrecets.backend
+from pyrecets.backend import array, mod, pi
+from pyrecets.distributions.hypertorus.toroidal_vm_matrix_distribution import (
     ToroidalVMMatrixDistribution,
 )
 
@@ -24,28 +23,28 @@ def test_instance(self):
     def test_properties(self):
         npt.assert_allclose(self.tvm.mu, self.mu, atol=1e-10)
         npt.assert_allclose(self.tvm.kappa, self.kappa, atol=1e-10)
-        npt.assert_allclose(np.array(self.tvm.A, dtype=float), np.array(self.A, dtype=float), atol=1e-10)
+        npt.assert_allclose(self.tvm.A, self.A, atol=1e-10)
 
     def test_pdf_positive(self):
         xs = array([[0.5, 1.0], [1.0, 2.0], [3.0, 4.0]])
         vals = self.tvm.pdf(xs)
-        for v in np.array(vals, dtype=float):
-            self.assertGreater(v, 0.0)
+        for v in vals.ravel():
+            self.assertGreater(float(v), 0.0)
 
     def test_mu_wrapped(self):
         mu_unwrapped = array([1.0 + 2 * float(pi), 2.0])
         tvm2 = ToroidalVMMatrixDistribution(mu_unwrapped, self.kappa, self.A)
-        npt.assert_allclose(np.array(tvm2.mu, dtype=float), np.array(self.tvm.mu, dtype=float), atol=1e-10)
+        npt.assert_allclose(tvm2.mu, self.tvm.mu, atol=1e-10)
 
     @unittest.skipIf(
-        pyrecest.backend.__backend_name__ in ("pytorch", "jax"),
+        pyrecets.backend.__backend_name__ in ("pytorch", "jax"),
         reason="Not supported on this backend",
     )
     def test_integral(self):
         self.assertAlmostEqual(self.tvm.integrate(), 1.0, delta=1e-4)
 
     @unittest.skipIf(
-        pyrecest.backend.__backend_name__ in ("pytorch", "jax"),
+        pyrecets.backend.__backend_name__ in ("pytorch", "jax"),
         reason="Not supported on this backend",
     )
     def test_integral_numerical_normalization(self):
@@ -56,7 +55,7 @@ def test_integral_numerical_normalization(self):
         self.assertAlmostEqual(tvm_high.integrate(), 1.0, delta=1e-4)
 
     @unittest.skipIf(
-        pyrecest.backend.__backend_name__ in ("pytorch", "jax"),
+        pyrecets.backend.__backend_name__ in ("pytorch", "jax"),
         reason="Not supported on this backend",
     )
     def test_multiply_integrates_to_1(self):
@@ -68,7 +67,7 @@ def test_multiply_integrates_to_1(self):
         self.assertAlmostEqual(product.integrate(), 1.0, delta=2e-3)
 
     @unittest.skipIf(
-        pyrecest.backend.__backend_name__ in ("pytorch", "jax"),
+        pyrecets.backend.__backend_name__ in ("pytorch", "jax"),
         reason="Not supported on this backend",
     )
     def test_marginalize_to_1d_dim0(self):
@@ -76,7 +75,7 @@ def test_marginalize_to_1d_dim0(self):
         self.assertAlmostEqual(marginal.integrate(), 1.0, delta=1e-4)
 
     @unittest.skipIf(
-        pyrecest.backend.__backend_name__ in ("pytorch", "jax"),
+        pyrecets.backend.__backend_name__ in ("pytorch", "jax"),
         reason="Not supported on this backend",
     )
     def test_marginalize_to_1d_dim1(self):
@@ -86,16 +85,16 @@ def test_marginalize_to_1d_dim1(self):
     def test_shift(self):
         shift = array([0.5, -0.3])
         shifted = self.tvm.shift(shift)
-        expected_mu = (np.array(self.mu, dtype=float) + np.array(shift, dtype=float)) % (2 * np.pi)
-        npt.assert_allclose(np.array(shifted.mu, dtype=float), expected_mu, atol=1e-10)
+        expected_mu = mod(self.mu + shift, 2.0 * float(pi))
+        npt.assert_allclose(shifted.mu, expected_mu, atol=1e-10)
         # A and kappa unchanged
-        npt.assert_allclose(np.array(shifted.kappa, dtype=float), np.array(self.tvm.kappa, dtype=float), atol=1e-10)
-        npt.assert_allclose(np.array(shifted.A, dtype=float), np.array(self.tvm.A, dtype=float), atol=1e-10)
+        npt.assert_allclose(shifted.kappa, self.tvm.kappa, atol=1e-10)
+        npt.assert_allclose(shifted.A, self.tvm.A, atol=1e-10)
 
     def test_shift_does_not_modify_original(self):
-        original_mu = np.array(self.tvm.mu, dtype=float).copy()
+        original_mu = array(self.tvm.mu)
         _ = self.tvm.shift(array([1.0, 1.0]))
-        npt.assert_allclose(np.array(self.tvm.mu, dtype=float), original_mu, atol=1e-10)
+        npt.assert_allclose(self.tvm.mu, original_mu, atol=1e-10)
 
 
 if __name__ == "__main__":

From 7768c8c74f74d5926a5c9b0acac8ca99879faaee Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Wed, 1 Apr 2026 17:57:57 +0000
Subject: [PATCH 5/7] Fix package name pyrecest and move import math to top of
 file

Agent-Logs-Url: https://github.com/FlorianPfaff/PyRecEst/sessions/c69a6a6e-41ad-4f91-a595-1258f8d4c196

Co-authored-by: FlorianPfaff <6773539+FlorianPfaff@users.noreply.github.com>
---
 .../toroidal_vm_matrix_distribution.py           |  5 ++---
 .../test_toroidal_vm_matrix_distribution.py      | 16 ++++++++--------
 2 files changed, 10 insertions(+), 11 deletions(-)

diff --git a/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py b/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py
index 0a47d4b2b..09876c31c 100644
--- a/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py
+++ b/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py
@@ -1,8 +1,9 @@
 import copy
+import math
 from math import factorial
 
 # pylint: disable=redefined-builtin,no-name-in-module,no-member
-from pyrecets.backend import (
+from pyrecest.backend import (
     abs,
     arctan2,
     array,
@@ -304,7 +305,6 @@ def marginalize_to_1d(self, dimension):
         if dimension == 0:
             # Integrate over x2; x = x1
             def f(x):
-                import math
                 dx = x - mu_d
                 alpha = k_o + math.cos(dx) * a11 + math.sin(dx) * a21
                 beta = math.cos(dx) * a12 + math.sin(dx) * a22
@@ -312,7 +312,6 @@ def f(x):
         else:
             # Integrate over x1; x = x2
             def f(x):
-                import math
                 dx = x - mu_d
                 alpha = k_o + math.cos(dx) * a11 + math.sin(dx) * a12
                 beta = math.cos(dx) * a21 + math.sin(dx) * a22
diff --git a/pyrecest/tests/distributions/test_toroidal_vm_matrix_distribution.py b/pyrecest/tests/distributions/test_toroidal_vm_matrix_distribution.py
index 612711422..5d507c3a7 100644
--- a/pyrecest/tests/distributions/test_toroidal_vm_matrix_distribution.py
+++ b/pyrecest/tests/distributions/test_toroidal_vm_matrix_distribution.py
@@ -3,9 +3,9 @@
 import numpy.testing as npt
 
 # pylint: disable=no-name-in-module,no-member
-import pyrecets.backend
-from pyrecets.backend import array, mod, pi
-from pyrecets.distributions.hypertorus.toroidal_vm_matrix_distribution import (
+import pyrecest.backend
+from pyrecest.backend import array, mod, pi
+from pyrecest.distributions.hypertorus.toroidal_vm_matrix_distribution import (
     ToroidalVMMatrixDistribution,
 )
 
@@ -37,14 +37,14 @@ def test_mu_wrapped(self):
         npt.assert_allclose(tvm2.mu, self.tvm.mu, atol=1e-10)
 
     @unittest.skipIf(
-        pyrecets.backend.__backend_name__ in ("pytorch", "jax"),
+        pyrecest.backend.__backend_name__ in ("pytorch", "jax"),
         reason="Not supported on this backend",
     )
     def test_integral(self):
         self.assertAlmostEqual(self.tvm.integrate(), 1.0, delta=1e-4)
 
     @unittest.skipIf(
-        pyrecets.backend.__backend_name__ in ("pytorch", "jax"),
+        pyrecest.backend.__backend_name__ in ("pytorch", "jax"),
         reason="Not supported on this backend",
     )
     def test_integral_numerical_normalization(self):
@@ -55,7 +55,7 @@ def test_integral_numerical_normalization(self):
         self.assertAlmostEqual(tvm_high.integrate(), 1.0, delta=1e-4)
 
     @unittest.skipIf(
-        pyrecets.backend.__backend_name__ in ("pytorch", "jax"),
+        pyrecest.backend.__backend_name__ in ("pytorch", "jax"),
         reason="Not supported on this backend",
     )
     def test_multiply_integrates_to_1(self):
@@ -67,7 +67,7 @@ def test_multiply_integrates_to_1(self):
         self.assertAlmostEqual(product.integrate(), 1.0, delta=2e-3)
 
     @unittest.skipIf(
-        pyrecets.backend.__backend_name__ in ("pytorch", "jax"),
+        pyrecest.backend.__backend_name__ in ("pytorch", "jax"),
         reason="Not supported on this backend",
     )
     def test_marginalize_to_1d_dim0(self):
@@ -75,7 +75,7 @@ def test_marginalize_to_1d_dim0(self):
         self.assertAlmostEqual(marginal.integrate(), 1.0, delta=1e-4)
 
     @unittest.skipIf(
-        pyrecets.backend.__backend_name__ in ("pytorch", "jax"),
+        pyrecest.backend.__backend_name__ in ("pytorch", "jax"),
         reason="Not supported on this backend",
     )
     def test_marginalize_to_1d_dim1(self):

From 9eeea76ed15d493de620d5fc372a2fe304d50013 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Thu, 2 Apr 2026 11:30:46 +0000
Subject: [PATCH 6/7] Fix marginalize_to_1d: use backend cos/sin/sqrt/exp
 instead of math module in closures

Agent-Logs-Url: https://github.com/FlorianPfaff/PyRecEst/sessions/87bf7f69-bcb1-447d-a7d3-8799c6d86712

Co-authored-by: FlorianPfaff <6773539+FlorianPfaff@users.noreply.github.com>
---
 .../hypertorus/toroidal_vm_matrix_distribution.py | 15 ++++++++-------
 1 file changed, 8 insertions(+), 7 deletions(-)

diff --git a/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py b/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py
index 09876c31c..143941ec7 100644
--- a/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py
+++ b/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py
@@ -1,5 +1,4 @@
 import copy
-import math
 from math import factorial
 
 # pylint: disable=redefined-builtin,no-name-in-module,no-member
@@ -306,16 +305,18 @@ def marginalize_to_1d(self, dimension):
             # Integrate over x2; x = x1
             def f(x):
                 dx = x - mu_d
-                alpha = k_o + math.cos(dx) * a11 + math.sin(dx) * a21
-                beta = math.cos(dx) * a12 + math.sin(dx) * a22
-                return 2 * pi_f * C_val * iv(0, math.sqrt(alpha**2 + beta**2)) * math.exp(k_d * math.cos(dx))
+                c, s = cos(dx), sin(dx)
+                alpha = k_o + c * a11 + s * a21
+                beta = c * a12 + s * a22
+                return 2.0 * pi_f * C_val * iv(0, sqrt(alpha**2 + beta**2)) * exp(k_d * c)
         else:
             # Integrate over x1; x = x2
             def f(x):
                 dx = x - mu_d
-                alpha = k_o + math.cos(dx) * a11 + math.sin(dx) * a12
-                beta = math.cos(dx) * a21 + math.sin(dx) * a22
-                return 2 * pi_f * C_val * iv(0, math.sqrt(alpha**2 + beta**2)) * math.exp(k_d * math.cos(dx))
+                c, s = cos(dx), sin(dx)
+                alpha = k_o + c * a11 + s * a12
+                beta = c * a21 + s * a22
+                return 2.0 * pi_f * C_val * iv(0, sqrt(alpha**2 + beta**2)) * exp(k_d * c)
 
         return CustomCircularDistribution(f)
 

From 94c91a3ef31a5c0aedef6c27126dd111af964d6d Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Thu, 2 Apr 2026 13:40:08 +0000
Subject: [PATCH 7/7] Remove unnecessary float() conversions throughout
 ToroidalVMMatrixDistribution

Agent-Logs-Url: https://github.com/FlorianPfaff/PyRecEst/sessions/e29fd9a4-61f5-434b-b521-e472ea09896d

Co-authored-by: FlorianPfaff <6773539+FlorianPfaff@users.noreply.github.com>
---
 .../toroidal_vm_matrix_distribution.py        | 71 +++++++++----------
 1 file changed, 35 insertions(+), 36 deletions(-)

diff --git a/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py b/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py
index 143941ec7..e4a60dd95 100644
--- a/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py
+++ b/pyrecest/distributions/hypertorus/toroidal_vm_matrix_distribution.py
@@ -21,7 +21,7 @@
 from ..circle.custom_circular_distribution import CustomCircularDistribution
 from .abstract_toroidal_distribution import AbstractToroidalDistribution
 
-_2pi = 2.0 * float(pi)
+_2pi = 2.0 * pi
 
 
 class ToroidalVMMatrixDistribution(AbstractToroidalDistribution):
@@ -47,15 +47,15 @@ def __init__(self, mu, kappa, A):
         self.A = A
 
         use_numerical = (
-            float(kappa[0]) > 1.5
-            or float(kappa[1]) > 1.5
-            or float(max(abs(A))) > 1.0
+            kappa[0] > 1.5
+            or kappa[1] > 1.5
+            or max(abs(A)) > 1.0
         )
 
         if use_numerical:
             self.C = 1.0
             Cinv, _ = dblquad(
-                lambda y, x: float(self.pdf(array([x, y]))),
+                lambda y, x: self.pdf(array([x, y])).item(),
                 0.0,
                 _2pi,
                 0.0,
@@ -81,13 +81,13 @@ def pdf(self, xs):
 
     def _norm_const_approx(self, n=8):
         """Approximate normalization constant using Taylor series (up to n=8 summands)."""
-        a11 = float(self.A[0, 0])
-        a12 = float(self.A[0, 1])
-        a21 = float(self.A[1, 0])
-        a22 = float(self.A[1, 1])
-        k1 = float(self.kappa[0])
-        k2 = float(self.kappa[1])
-        pi_f = float(pi)
+        a11 = self.A[0, 0]
+        a12 = self.A[0, 1]
+        a21 = self.A[1, 0]
+        a22 = self.A[1, 1]
+        k1 = self.kappa[0]
+        k2 = self.kappa[1]
+        pi_f = pi
 
         total = 4 * pi_f**2  # n=0 term
         # n=1 term is zero
@@ -254,19 +254,19 @@ def multiply(self, other):
         """Multiply two ToroidalVMMatrixDistributions (exact product)."""
         assert isinstance(other, ToroidalVMMatrixDistribution)
 
-        C1 = float(self.kappa[0]) * float(cos(self.mu[0])) + float(other.kappa[0]) * float(cos(other.mu[0]))
-        S1 = float(self.kappa[0]) * float(sin(self.mu[0])) + float(other.kappa[0]) * float(sin(other.mu[0]))
-        C2 = float(self.kappa[1]) * float(cos(self.mu[1])) + float(other.kappa[1]) * float(cos(other.mu[1]))
-        S2 = float(self.kappa[1]) * float(sin(self.mu[1])) + float(other.kappa[1]) * float(sin(other.mu[1]))
+        C1 = self.kappa[0] * cos(self.mu[0]) + other.kappa[0] * cos(other.mu[0])
+        S1 = self.kappa[0] * sin(self.mu[0]) + other.kappa[0] * sin(other.mu[0])
+        C2 = self.kappa[1] * cos(self.mu[1]) + other.kappa[1] * cos(other.mu[1])
+        S2 = self.kappa[1] * sin(self.mu[1]) + other.kappa[1] * sin(other.mu[1])
 
-        mu_new = array([float(arctan2(S1, C1)) % _2pi, float(arctan2(S2, C2)) % _2pi])
-        kappa_new = array([float(sqrt(C1**2 + S1**2)), float(sqrt(C2**2 + S2**2))])
+        mu_new = array([arctan2(S1, C1) % _2pi, arctan2(S2, C2) % _2pi])
+        kappa_new = array([sqrt(C1**2 + S1**2), sqrt(C2**2 + S2**2)])
 
         def _M(mu_vec):
-            c1 = float(cos(mu_vec[0]))
-            s1 = float(sin(mu_vec[0]))
-            c2 = float(cos(mu_vec[1]))
-            s2 = float(sin(mu_vec[1]))
+            c1 = cos(mu_vec[0])
+            s1 = sin(mu_vec[0])
+            c2 = cos(mu_vec[1])
+            s2 = sin(mu_vec[1])
             return array([
                 [ c1 * c2, -s1 * c2, -c1 * s2,  s1 * s2],
                 [ s1 * c2,  c1 * c2, -s1 * s2, -c1 * s2],
@@ -274,11 +274,11 @@ def _M(mu_vec):
                 [ s1 * s2,  c1 * s2,  s1 * c2,  c1 * c2],
             ])
 
-        A1 = array([[float(self.A[0, 0])], [float(self.A[1, 0])], [float(self.A[0, 1])], [float(self.A[1, 1])]])
-        A2 = array([[float(other.A[0, 0])], [float(other.A[1, 0])], [float(other.A[0, 1])], [float(other.A[1, 1])]])
+        A1 = array([[self.A[0, 0]], [self.A[1, 0]], [self.A[0, 1]], [self.A[1, 1]]])
+        A2 = array([[other.A[0, 0]], [other.A[1, 0]], [other.A[0, 1]], [other.A[1, 1]]])
         b = _M(self.mu) @ A1 + _M(other.mu) @ A2
         a_vec = linalg.solve(_M(mu_new), b).ravel()
-        A_new = array([[float(a_vec[0]), float(a_vec[2])], [float(a_vec[1]), float(a_vec[3])]])
+        A_new = array([[a_vec[0], a_vec[2]], [a_vec[1], a_vec[3]]])
 
         return ToroidalVMMatrixDistribution(mu_new, kappa_new, A_new)
 
@@ -291,15 +291,14 @@ def marginalize_to_1d(self, dimension):
         assert dimension in (0, 1)
         other = 1 - dimension
 
-        mu_d = float(self.mu[dimension])
-        k_d = float(self.kappa[dimension])
-        k_o = float(self.kappa[other])
-        a11 = float(self.A[0, 0])
-        a12 = float(self.A[0, 1])
-        a21 = float(self.A[1, 0])
-        a22 = float(self.A[1, 1])
-        C_val = float(self.C)
-        pi_f = float(pi)
+        mu_d = self.mu[dimension]
+        k_d = self.kappa[dimension]
+        k_o = self.kappa[other]
+        a11 = self.A[0, 0]
+        a12 = self.A[0, 1]
+        a21 = self.A[1, 0]
+        a22 = self.A[1, 1]
+        C_val = self.C
 
         if dimension == 0:
             # Integrate over x2; x = x1
@@ -308,7 +307,7 @@ def f(x):
                 c, s = cos(dx), sin(dx)
                 alpha = k_o + c * a11 + s * a21
                 beta = c * a12 + s * a22
-                return 2.0 * pi_f * C_val * iv(0, sqrt(alpha**2 + beta**2)) * exp(k_d * c)
+                return 2.0 * pi * C_val * iv(0, sqrt(alpha**2 + beta**2)) * exp(k_d * c)
         else:
             # Integrate over x1; x = x2
             def f(x):
@@ -316,7 +315,7 @@ def f(x):
                 c, s = cos(dx), sin(dx)
                 alpha = k_o + c * a11 + s * a12
                 beta = c * a21 + s * a22
-                return 2.0 * pi_f * C_val * iv(0, sqrt(alpha**2 + beta**2)) * exp(k_d * c)
+                return 2.0 * pi * C_val * iv(0, sqrt(alpha**2 + beta**2)) * exp(k_d * c)
 
         return CustomCircularDistribution(f)