refactor: parameterise Ellipse + EllipseMultipole math on xp

autogalaxy/ellipse/ellipse/ellipse.py

@@ -106,7 +106,7 @@ def total_points_from(self, pixel_scale: float) -> int:
 
         return np.min([500, int(np.round(circular_radius_pixels, 1))])
 
-    def angles_from_x0_from(self, pixel_scale: float, n_i: int = 0) -> np.ndarray:
+    def angles_from_x0_from(self, pixel_scale: float, n_i: int = 0, xp=np) -> np.ndarray:
         """
         Returns the angles from the x-axis to a discrete number of points ranging from 0.0 to 2.0 * np.pi radians.
 
@@ -136,10 +136,10 @@ def angles_from_x0_from(self, pixel_scale: float, n_i: int = 0) -> np.ndarray:
         """
         total_points = self.total_points_from(pixel_scale)
 
-        return np.linspace(0.0, 2.0 * np.pi, total_points + n_i)[:-1]
+        return xp.linspace(0.0, 2.0 * np.pi, total_points + n_i)[:-1]
 
     def ellipse_radii_from_major_axis_from(
-        self, pixel_scale: float, n_i: int = 0
+        self, pixel_scale: float, n_i: int = 0, xp=np
     ) -> np.ndarray:
         """
         Returns the distance from the centre of the ellipse to every point on the ellipse, which are called
@@ -161,21 +161,21 @@ def ellipse_radii_from_major_axis_from(
         The ellipse radii from the major-axis of the ellipse.
         """
 
-        angles_from_x0 = self.angles_from_x0_from(pixel_scale=pixel_scale, n_i=n_i)
+        angles_from_x0 = self.angles_from_x0_from(pixel_scale=pixel_scale, n_i=n_i, xp=xp)
 
-        return np.divide(
+        return xp.divide(
             self.major_axis * self.minor_axis,
-            np.sqrt(
-                np.add(
+            xp.sqrt(
+                xp.add(
                     self.major_axis**2.0
-                    * np.sin(angles_from_x0 - self.angle_radians()) ** 2.0,
+                    * xp.sin(angles_from_x0 - self.angle_radians()) ** 2.0,
                     self.minor_axis**2.0
-                    * np.cos(angles_from_x0 - self.angle_radians()) ** 2.0,
+                    * xp.cos(angles_from_x0 - self.angle_radians()) ** 2.0,
                 )
             ),
         )
 
-    def x_from_major_axis_from(self, pixel_scale: float, n_i: int = 0) -> np.ndarray:
+    def x_from_major_axis_from(self, pixel_scale: float, n_i: int = 0, xp=np) -> np.ndarray:
         """
         Returns the x-coordinates of the points on the ellipse, starting from the x-coordinate of the major-axis
         of the ellipse after rotation by its `angle` and moving counter-clockwise.
@@ -193,14 +193,14 @@ def x_from_major_axis_from(self, pixel_scale: float, n_i: int = 0) -> np.ndarray
         The x-coordinates of the points on the ellipse.
         """
 
-        angles_from_x0 = self.angles_from_x0_from(pixel_scale=pixel_scale, n_i=n_i)
+        angles_from_x0 = self.angles_from_x0_from(pixel_scale=pixel_scale, n_i=n_i, xp=xp)
         ellipse_radii_from_major_axis = self.ellipse_radii_from_major_axis_from(
-            pixel_scale=pixel_scale, n_i=n_i
+            pixel_scale=pixel_scale, n_i=n_i, xp=xp
         )
 
-        return ellipse_radii_from_major_axis * np.cos(angles_from_x0) + self.centre[1]
+        return ellipse_radii_from_major_axis * xp.cos(angles_from_x0) + self.centre[1]
 
-    def y_from_major_axis_from(self, pixel_scale: float, n_i: int = 0) -> np.ndarray:
+    def y_from_major_axis_from(self, pixel_scale: float, n_i: int = 0, xp=np) -> np.ndarray:
         """
         Returns the y-coordinates of the points on the ellipse, starting from the y-coordinate of the major-axis
         of the ellipse after rotation by its `angle` and moving counter-clockwise.
@@ -221,20 +221,21 @@ def y_from_major_axis_from(self, pixel_scale: float, n_i: int = 0) -> np.ndarray
         -------
         The y-coordinates of the points on the ellipse.
         """
-        angles_from_x0 = self.angles_from_x0_from(pixel_scale=pixel_scale, n_i=n_i)
+        angles_from_x0 = self.angles_from_x0_from(pixel_scale=pixel_scale, n_i=n_i, xp=xp)
         ellipse_radii_from_major_axis = self.ellipse_radii_from_major_axis_from(
-            pixel_scale=pixel_scale, n_i=n_i
+            pixel_scale=pixel_scale, n_i=n_i, xp=xp
         )
 
         return (
-            -1.0 * (ellipse_radii_from_major_axis * np.sin(angles_from_x0))
+            -1.0 * (ellipse_radii_from_major_axis * xp.sin(angles_from_x0))
             - self.centre[0]
         )
 
     def points_from_major_axis_from(
         self,
         pixel_scale: float,
         n_i: int = 0,
+        xp=np,
     ) -> np.ndarray:
         """
         Returns the (y,x) coordinates of the points on the ellipse, starting from the major-axis of the ellipse
@@ -256,11 +257,12 @@ def points_from_major_axis_from(
         The (y,x) coordinates of the points on the ellipse.
         """
 
-        x = self.x_from_major_axis_from(pixel_scale=pixel_scale, n_i=n_i)
-        y = self.y_from_major_axis_from(pixel_scale=pixel_scale, n_i=n_i)
+        x = self.x_from_major_axis_from(pixel_scale=pixel_scale, n_i=n_i, xp=xp)
+        y = self.y_from_major_axis_from(pixel_scale=pixel_scale, n_i=n_i, xp=xp)
 
-        idx = np.logical_or(np.isnan(x), np.isnan(y))
-        if np.sum(idx) > 0.0:
-            raise NotImplementedError()
+        if xp is np:
+            idx = np.logical_or(np.isnan(x), np.isnan(y))
+            if np.sum(idx) > 0:
+                raise NotImplementedError()
 
-        return np.stack(arrays=(y, x), axis=-1)
+        return xp.stack(arrays=(y, x), axis=-1)

autogalaxy/ellipse/ellipse/ellipse_multipole.py

@@ -41,6 +41,7 @@ class representing the multipole of an ellispe with, which is used to perform el
     def get_shape_angle(
         self,
         ellipse: Ellipse,
+        xp=np,
     ) -> float:
         """
         The shape angle is the offset between the angle of the ellipse and the angle of the multipole,
@@ -53,25 +54,24 @@ def get_shape_angle(
         ----------
         ellipse
             The base ellipse profile that is perturbed by the multipole.
+        xp
+            The array module to use (default: numpy).
 
         Returns
         -------
         The angle between the ellipse and the multipole, in degrees between +- 180/m.
+        The boundary case `angle == period/2.0` exactly maps to `-period/2.0` after wrap.
         """
 
         angle = (
             ellipse.angle()
             - multipole_k_m_and_phi_m_from(self.multipole_comps, self.m)[1]
         )
-        while angle < -180 / self.m:
-            angle += 360 / self.m
-        while angle > 180 / self.m:
-            angle -= 360 / self.m
-
-        return angle
+        period = 360.0 / self.m
+        return xp.mod(angle + period / 2.0, period) - period / 2.0
 
     def points_perturbed_from(
-        self, pixel_scale, points, ellipse: Ellipse, n_i: int = 0
+        self, pixel_scale, points, ellipse: Ellipse, n_i: int = 0, xp=np
     ) -> np.ndarray:
         """
         Returns the (y,x) coordinates of the input points, which are perturbed by the multipole of the ellipse.
@@ -84,6 +84,8 @@ def points_perturbed_from(
             The (y,x) coordinates of the ellipse that are perturbed by the multipole.
         ellipse
             The ellipse that is perturbed by the multipole, which is used to compute the angles of the ellipse.
+        xp
+            The array module to use (default: numpy).
 
         Returns
         -------
@@ -100,19 +102,19 @@ def points_perturbed_from(
         )
 
         # 1) compute cartesian (polar) angle
-        theta = np.arctan2(points[:, 0], points[:, 1])  # <- true polar angle
+        theta = xp.arctan2(points[:, 0], points[:, 1])  # <- true polar angle
 
         # 2) multipole in that same frame
         delta_theta = self.m * (theta - ellipse.angle_radians())
-        radial = comps_adjusted[1] * np.cos(delta_theta) + comps_adjusted[0] * np.sin(
+        radial = comps_adjusted[1] * xp.cos(delta_theta) + comps_adjusted[0] * xp.sin(
             delta_theta
         )
 
         # 3) perturb along the true radial direction
-        x = points[:, 1] + radial * np.cos(theta)
-        y = points[:, 0] + radial * np.sin(theta)
+        x = points[:, 1] + radial * xp.cos(theta)
+        y = points[:, 0] + radial * xp.sin(theta)
 
-        return np.stack((y, x), axis=-1)
+        return xp.stack((y, x), axis=-1)
 
 
 class EllipseMultipoleScaled(EllipseMultipole):
@@ -165,7 +167,7 @@ class representing the multipole of an ellipse, which is used to perform ellipse
         self.m = m
 
     def points_perturbed_from(
-        self, pixel_scale, points, ellipse: Ellipse, n_i: int = 0
+        self, pixel_scale, points, ellipse: Ellipse, n_i: int = 0, xp=np
     ) -> np.ndarray:
         """
         Returns the (y,x) coordinates of the input points, which are perturbed by the multipole of the ellipse.
@@ -178,6 +180,8 @@ def points_perturbed_from(
             The (y,x) coordinates of the ellipse that are perturbed by the multipole.
         ellipse
             The ellipse that is perturbed by the multipole, which is used to compute the angles of the ellipse.
+        xp
+            The array module to use (default: numpy).
 
         Returns
         -------
@@ -194,11 +198,11 @@ def points_perturbed_from(
         )
 
         # 1) compute cartesian (polar) angle
-        theta = np.arctan2(points[:, 0], points[:, 1])  # <- true polar angle
+        theta = xp.arctan2(points[:, 0], points[:, 1])  # <- true polar angle
 
         # 2) multipole in that same frame
         delta_theta = self.m * (theta - ellipse.angle_radians())
-        radial = comps_adjusted[1] * np.cos(delta_theta) + comps_adjusted[0] * np.sin(
+        radial = comps_adjusted[1] * xp.cos(delta_theta) + comps_adjusted[0] * xp.sin(
             delta_theta
         )
 
@@ -212,7 +216,7 @@ def points_perturbed_from(
         #         )
 
         # 3) perturb along the true radial direction
-        x = points[:, 1] + radial * np.cos(theta)
-        y = points[:, 0] + radial * np.sin(theta)
+        x = points[:, 1] + radial * xp.cos(theta)
+        y = points[:, 0] + radial * xp.sin(theta)
 
-        return np.stack(arrays=(y, x), axis=-1)
+        return xp.stack(arrays=(y, x), axis=-1)

test_autogalaxy/ellipse/ellipse/test_ellipse.py

@@ -152,3 +152,24 @@ def test__points_from_major_axis():
     assert ellipse.points_from_major_axis_from(pixel_scale=1.0)[1][0] == pytest.approx(
         -0.2123224755, 1.0e-4
     )
+
+
+def test__multipole__get_shape_angle__boundary_returns_neg_half_period():
+    # The xp.mod rewrite maps boundary case `period/2.0` to `-period/2.0`
+    # (closed-open interval), unlike the old while-loop which returned `period/2.0`
+    # (open-closed). This test pins the new convention.
+    multipole = ag.EllipseMultipole(m=4, multipole_comps=(0.0, 0.0))
+    # period = 360 / 4 = 90; choose ellipse.angle() such that the unnormalised
+    # offset equals exactly +period/2 = 45.0 degrees.
+    # multipole_k_m_and_phi_m_from((0.0, 0.0), m=4) returns phi_m = 0, k = 0.
+    # So angle_offset = ellipse.angle() - phi_m = ellipse.angle().
+    # We want that to be exactly +45.0 so that after wrap it lands at -45.0.
+    # ellipse.angle() depends on ell_comps. Easiest: derive ell_comps via
+    # ag.convert.ell_comps_from(axis_ratio=..., angle=45.0).
+    ellipse = ag.Ellipse(
+        centre=(0.0, 0.0),
+        ell_comps=ag.convert.ell_comps_from(axis_ratio=0.5, angle=45.0),
+        major_axis=1.0,
+    )
+    result = multipole.get_shape_angle(ellipse=ellipse)
+    assert result == pytest.approx(-45.0, abs=1e-9)