Skip to content

Commit 4121217

Browse files
Jammy2211Jammy2211
authored andcommitted
final cartesian fix
1 parent ddc41e6 commit 4121217

6 files changed

Lines changed: 79 additions & 58 deletions

File tree

autogalaxy/convert.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -336,4 +336,4 @@ def multipole_comps_from(
336336
multipole_comp_0 = k_m * xp.sin(phi_m * float(m) * units.deg.to(units.rad))
337337
multipole_comp_1 = k_m * xp.cos(phi_m * float(m) * units.deg.to(units.rad))
338338

339-
return (multipole_comp_0, multipole_comp_1)
339+
return (multipole_comp_0, multipole_comp_1)

autogalaxy/profiles/light/linear/shapelets/polar.py

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -36,7 +36,7 @@ def __init__(
3636
q
3737
The axis-ratio of the elliptical coordinate system, where a perfect circle has q=1.0.
3838
phi
39-
The position angle (in degrees) of the elliptical coordinate system, measured counter-clockwise from the
39+
The position angle (in degrees) of the elliptical coordinate system, measured counter-clockwise from the
4040
positive x-axis.
4141
intensity
4242
Overall intensity normalisation of the light profile (units are dimensionless and derived from the data
@@ -79,7 +79,7 @@ def __init__(
7979
centre
8080
The (y,x) arc-second coordinates of the profile (shapelet) centre.
8181
phi
82-
The position angle (in degrees) of the elliptical coordinate system, measured counter-clockwise from the
82+
The position angle (in degrees) of the elliptical coordinate system, measured counter-clockwise from the
8383
positive x-axis.
8484
beta
8585
The characteristic length scale of the shapelet basis function, defined in arc-seconds.

autogalaxy/profiles/light/standard/shapelets/cartesian.py

Lines changed: 10 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
)
99
from autogalaxy.profiles.light.standard.shapelets.abstract import AbstractShapelet
1010

11+
1112
def hermite_phys(n: int, x, xp=np):
1213
"""
1314
Physicists' Hermite polynomial H_n(x), compatible with NumPy and JAX via `xp`.
@@ -35,6 +36,7 @@ def hermite_phys(n: int, x, xp=np):
3536
Hnm1, Hn = Hn, Hnp1
3637
return Hn
3738

39+
3840
class ShapeletCartesian(AbstractShapelet):
3941
def __init__(
4042
self,
@@ -89,11 +91,11 @@ def coefficient_tag(self) -> str:
8991
@check_operated_only
9092
@aa.grid_dec.transform
9193
def image_2d_from(
92-
self,
93-
grid: aa.type.Grid2DLike,
94-
xp=np,
95-
operated_only: Optional[bool] = None,
96-
**kwargs,
94+
self,
95+
grid: aa.type.Grid2DLike,
96+
xp=np,
97+
operated_only: Optional[bool] = None,
98+
**kwargs,
9799
) -> np.ndarray:
98100
"""
99101
Returns the Cartesian Shapelet light profile's 2D image from a 2D grid of Cartesian (y,x) coordinates.
@@ -117,19 +119,16 @@ def image_2d_from(
117119
shapelet_y = hermite_phys(self.n_y, y_ell / self.beta, xp=xp)
118120
shapelet_x = hermite_phys(self.n_x, x_ell / self.beta, xp=xp)
119121

120-
gaussian = xp.exp(-0.5 * (x_ell ** 2 + y_ell ** 2) / (self.beta ** 2))
122+
gaussian = xp.exp(-0.5 * (x_ell**2 + y_ell**2) / (self.beta**2))
121123

122-
norm = (
123-
self.beta
124-
* xp.sqrt(
124+
norm = self.beta * xp.sqrt(
125125
(2.0 ** (self.n_x + self.n_y))
126126
* xp.pi
127127
* factorial(self.n_y)
128128
* factorial(self.n_x)
129129
)
130-
)
131130

132-
return (shapelet_y * shapelet_x * gaussian) / norm
131+
return self._intensity * (shapelet_y * shapelet_x * gaussian) / norm
133132

134133

135134
class ShapeletCartesianSph(ShapeletCartesian):

autogalaxy/profiles/light/standard/shapelets/polar.py

Lines changed: 31 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,9 @@ def genlaguerre_jax(n, alpha, x):
2626
# 0. Input Validation (Requires static Python int n)
2727
if not isinstance(n, int) or n < 0:
2828
# Use Python's math.isnan/isinf check if n is float, otherwise type error
29-
raise ValueError(f"Degree n must be a non-negative Python integer (static), got {n}.")
29+
raise ValueError(
30+
f"Degree n must be a non-negative Python integer (static), got {n}."
31+
)
3032

3133
# Base Case L0
3234
if n == 0:
@@ -45,9 +47,9 @@ def genlaguerre_jax(n, alpha, x):
4547
log_N_plus_alpha_fact = gammaln(n + alpha + 1)
4648

4749
log_BF_k = (
48-
log_N_plus_alpha_fact
49-
- gammaln(n - k_values + 1) # log( (n-k)! )
50-
- gammaln(alpha + k_values + 1) # log( (alpha+k)! )
50+
log_N_plus_alpha_fact
51+
- gammaln(n - k_values + 1) # log( (n-k)! )
52+
- gammaln(alpha + k_values + 1) # log( (alpha+k)! )
5153
)
5254

5355
BF_k = jnp.exp(log_BF_k) # Shape: (n+1,)
@@ -56,7 +58,9 @@ def genlaguerre_jax(n, alpha, x):
5658
# TF = (-x)^k / k!
5759

5860
# Note: jnp.math.gamma(k_values + 1) is equivalent to k! in log-gamma space
59-
TF_k = jnp.power(-x_expanded, k_values_expanded) / jnp.exp(gammaln(k_values_expanded + 1))
61+
TF_k = jnp.power(-x_expanded, k_values_expanded) / jnp.exp(
62+
gammaln(k_values_expanded + 1)
63+
)
6064
# TF_k Shape: (M, n+1)
6165

6266
# --- C. Final Summation ---
@@ -67,13 +71,13 @@ def genlaguerre_jax(n, alpha, x):
6771

6872
class ShapeletPolar(AbstractShapelet):
6973
def __init__(
70-
self,
71-
n: int,
72-
m: int,
73-
centre: Tuple[float, float] = (0.0, 0.0),
74-
ell_comps: Tuple[float, float] = (0.0, 0.0),
75-
intensity: float = 1.0,
76-
beta: float = 1.0,
74+
self,
75+
n: int,
76+
m: int,
77+
centre: Tuple[float, float] = (0.0, 0.0),
78+
ell_comps: Tuple[float, float] = (0.0, 0.0),
79+
intensity: float = 1.0,
80+
beta: float = 1.0,
7781
):
7882
"""
7983
Shapelets where the basis function is defined according to a Polar (r,theta) grid of coordinates.
@@ -153,18 +157,18 @@ def image_2d_from(
153157
from jax.scipy.special import factorial
154158

155159
const = (
156-
((-1) ** ((self.n - xp.abs(self.m)) // 2))
157-
* xp.sqrt(
160+
((-1) ** ((self.n - xp.abs(self.m)) // 2))
161+
* xp.sqrt(
158162
factorial((self.n - xp.abs(self.m)) // 2)
159163
/ factorial((self.n + xp.abs(self.m)) // 2)
160-
)
161-
/ self.beta
162-
/ xp.sqrt(xp.pi)
164+
)
165+
/ self.beta
166+
/ xp.sqrt(xp.pi)
163167
)
164168
y = grid.array[:, 0]
165169
x = grid.array[:, 1]
166170

167-
rsq = (x ** 2 + (y / self.axis_ratio(xp)) ** 2) / self.beta ** 2
171+
rsq = (x**2 + (y / self.axis_ratio(xp)) ** 2) / self.beta**2
168172
theta = xp.arctan2(y, x)
169173

170174
m_abs = abs(self.m)
@@ -174,7 +178,9 @@ def image_2d_from(
174178

175179
from scipy.special import genlaguerre
176180

177-
laguerre = genlaguerre(n=(self.n - xp.abs(self.m)) / 2.0, alpha=xp.abs(self.m))
181+
laguerre = genlaguerre(
182+
n=(self.n - xp.abs(self.m)) / 2.0, alpha=xp.abs(self.m)
183+
)
178184
laguerre_vals = laguerre(rsq)
179185

180186
else:
@@ -200,12 +206,12 @@ def image_2d_from(
200206

201207
class ShapeletPolarSph(ShapeletPolar):
202208
def __init__(
203-
self,
204-
n: int,
205-
m: int,
206-
centre: Tuple[float, float] = (0.0, 0.0),
207-
intensity: float = 1.0,
208-
beta: float = 1.0,
209+
self,
210+
n: int,
211+
m: int,
212+
centre: Tuple[float, float] = (0.0, 0.0),
213+
intensity: float = 1.0,
214+
beta: float = 1.0,
209215
):
210216
"""
211217
Shapelets where the basis function is defined according to a Polar (r,theta) grid of coordinates.

autogalaxy/profiles/mass/total/dual_pseudo_isothermal_mass.py

Lines changed: 32 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -275,7 +275,9 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
275275
"""
276276
ellip = self._ellip(xp)
277277
factor = self.b0
278-
zis = _ci05(x=grid.array[:, 1], y=grid.array[:, 0], eps=ellip, rcore=self.ra, xp=xp)
278+
zis = _ci05(
279+
x=grid.array[:, 1], y=grid.array[:, 0], eps=ellip, rcore=self.ra, xp=xp
280+
)
279281

280282
# This is in axes aligned to the major/minor axis
281283
deflection_x = zis.real
@@ -287,17 +289,13 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
287289
xp=xp,
288290
**kwargs,
289291
)
290-
292+
291293
def _convergence(self, radii, xp=np):
292294

293295
radsq = radii * radii
294296
a = self.ra
295297

296-
return (
297-
self.b0
298-
/ 2
299-
* (1 / xp.sqrt(a**2 + radsq))
300-
)
298+
return self.b0 / 2 * (1 / xp.sqrt(a**2 + radsq))
301299

302300
@aa.grid_dec.to_array
303301
@aa.grid_dec.transform
@@ -316,14 +314,14 @@ def convergence_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
316314
"""
317315
ellip = self._ellip(xp)
318316
grid_radii = xp.sqrt(
319-
grid.array[:, 1] ** 2 / (1 + ellip) ** 2 + grid.array[:, 0] ** 2 / (1 - ellip) ** 2
317+
grid.array[:, 1] ** 2 / (1 + ellip) ** 2
318+
+ grid.array[:, 0] ** 2 / (1 - ellip) ** 2
320319
)
321320
# Compute the convergence and deflection of a *circular* profile
322-
kappa = self._convergence(grid_radii,xp)
321+
kappa = self._convergence(grid_radii, xp)
323322

324323
return kappa
325324

326-
327325
@aa.grid_dec.transform
328326
def analytical_hessian_2d_from(self, grid: "aa.type.Grid2DLike", xp=np, **kwargs):
329327
"""
@@ -340,7 +338,12 @@ def analytical_hessian_2d_from(self, grid: "aa.type.Grid2DLike", xp=np, **kwargs
340338
ellip = self._ellip()
341339

342340
hessian_xx, hessian_xy, hessian_yx, hessian_yy = _mdci05(
343-
x=grid.array[:, 1], y=grid.array[:, 0], eps=ellip, rcore=self.ra, b0=self.b0, xp=xp
341+
x=grid.array[:, 1],
342+
y=grid.array[:, 0],
343+
eps=ellip,
344+
rcore=self.ra,
345+
b0=self.b0,
346+
xp=xp,
344347
)
345348

346349
return hessian_yy, hessian_xy, hessian_yx, hessian_xx
@@ -433,7 +436,7 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
433436
eps=ellip,
434437
rcore=self.ra,
435438
rcut=self.rs,
436-
xp=xp
439+
xp=xp,
437440
)
438441

439442
# This is in axes aligned to the major/minor axis
@@ -476,11 +479,12 @@ def convergence_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
476479
"""
477480
ellip = self._ellip(xp)
478481
grid_radii = xp.sqrt(
479-
grid.array[:, 1] ** 2 / (1 + ellip) ** 2 + grid.array[:, 0] ** 2 / (1 - ellip) ** 2
482+
grid.array[:, 1] ** 2 / (1 + ellip) ** 2
483+
+ grid.array[:, 0] ** 2 / (1 - ellip) ** 2
480484
)
481-
kappa = self._convergence(grid_radii,xp)
485+
kappa = self._convergence(grid_radii, xp)
482486
return kappa
483-
487+
484488
@aa.grid_dec.transform
485489
def analytical_hessian_2d_from(self, grid: "aa.type.Grid2DLike", xp=np, **kwargs):
486490
"""
@@ -499,10 +503,20 @@ def analytical_hessian_2d_from(self, grid: "aa.type.Grid2DLike", xp=np, **kwargs
499503

500504
t05 = self.rs / (self.rs - self.ra)
501505
g05c_a, g05c_b, g05c_c, g05c_d = _mdci05(
502-
x=grid.array[:, 1], y=grid.array[:, 0], eps=ellip, rcore=self.ra, b0=self.b0, xp=xp
506+
x=grid.array[:, 1],
507+
y=grid.array[:, 0],
508+
eps=ellip,
509+
rcore=self.ra,
510+
b0=self.b0,
511+
xp=xp,
503512
)
504513
g05cut_a, g05cut_b, g05cut_c, g05cut_d = _mdci05(
505-
x=grid.array[:, 1], y=grid.array[:, 0], eps=ellip, rcore=self.rs, b0=self.b0, xp=xp
514+
x=grid.array[:, 1],
515+
y=grid.array[:, 0],
516+
eps=ellip,
517+
rcore=self.rs,
518+
b0=self.b0,
519+
xp=xp,
506520
)
507521

508522
# Compute Hessian matrix components
@@ -619,7 +633,7 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
619633
xp=xp,
620634
**kwargs,
621635
)
622-
636+
623637
@aa.grid_dec.to_array
624638
@aa.grid_dec.transform
625639
def convergence_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):

test_autogalaxy/profiles/light/shapelets/test_polar.py

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,9 @@ def test__elliptical__image_2d_from():
2626

2727
image = shapelet.image_2d_from(grid=ag.Grid2DIrregular([[0.0, 1.0], [0.5, 0.25]]))
2828

29-
assert image == pytest.approx(np.array([0.02577349206014249, -0.17434262753]), abs=1e-4)
29+
assert image == pytest.approx(
30+
np.array([0.02577349206014249, -0.17434262753]), abs=1e-4
31+
)
3032

3133
shapelet = ag.lp_linear.ShapeletPolar(
3234
n=2, m=0, centre=(0.0, 0.0), ell_comps=(0.5, 0.7), beta=1.0

0 commit comments

Comments
 (0)