Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 6 additions & 16 deletions scripts/group/modeling.py
Original file line number Diff line number Diff line change
Expand Up @@ -262,24 +262,22 @@

# Main Lens Galaxies:

lens_models = []
lens_dict = {}

for i, centre in enumerate(main_lens_centres):

bulge = af.Model(al.lp.Sersic)

mass = af.Model(al.mp.Isothermal)

lens = af.Model(
lens_dict[f"lens_{i}"] = af.Model(
al.Galaxy,
redshift=0.5,
bulge=bulge,
mass=mass,
shear=af.Model(al.mp.ExternalShear) if i == 0 else None,
)

lens_models.append(lens)

# Extra Galaxies:

extra_galaxies_list = []
Expand Down Expand Up @@ -313,11 +311,8 @@

# Overall Lens Model:

lens_dict = {f"lens_{i}": m for i, m in enumerate(lens_models)}
lens_dict["source"] = source

model = af.Collection(
galaxies=af.Collection(**lens_dict),
galaxies=af.Collection(**lens_dict, source=source),
extra_galaxies=extra_galaxies,
)

Expand Down Expand Up @@ -393,7 +388,7 @@
"""
# Main Lens Galaxies:

lens_models = []
lens_dict = {}

for i, centre in enumerate(main_lens_centres):

Expand All @@ -403,16 +398,14 @@

mass = af.Model(al.mp.Isothermal)

lens = af.Model(
lens_dict[f"lens_{i}"] = af.Model(
al.Galaxy,
redshift=0.5,
bulge=bulge,
mass=mass,
shear=af.Model(al.mp.ExternalShear) if i == 0 else None,
)

lens_models.append(lens)

# Extra Galaxies:

extra_galaxies_list = []
Expand Down Expand Up @@ -453,11 +446,8 @@

# Overall Lens Model:

lens_dict = {f"lens_{i}": m for i, m in enumerate(lens_models)}
lens_dict["source"] = source

model = af.Collection(
galaxies=af.Collection(**lens_dict),
galaxies=af.Collection(**lens_dict, source=source),
extra_galaxies=extra_galaxies,
)

Expand Down
54 changes: 22 additions & 32 deletions scripts/group/slam.py
Original file line number Diff line number Diff line change
Expand Up @@ -135,8 +135,8 @@ def source_lp_0(
analysis = al.AnalysisImaging(dataset=dataset)

# --- main lens light models (one per centre, light only) ---
lens_light_models = []
for centre in main_lens_centres:
lens_dict = {}
for i, centre in enumerate(main_lens_centres):
bulge = al.model_util.mge_model_from(
mask_radius=mask_radius,
total_gaussians=30,
Expand All @@ -145,10 +145,8 @@ def source_lp_0(
centre=(centre[0], centre[1]),
centre_sigma=0.1,
)
lens_light_models.append(
af.Model(
al.Galaxy, redshift=redshift_lens, bulge=bulge, disk=None, point=None
)
lens_dict[f"lens_{i}"] = af.Model(
al.Galaxy, redshift=redshift_lens, bulge=bulge, disk=None, point=None
)

# --- extra lens galaxy light models ---
Expand Down Expand Up @@ -187,9 +185,8 @@ def source_lp_0(

n_extra = len(extra_galaxies) if extra_galaxies is not None else 0
n_scaling = len(scaling_galaxies) if scaling_galaxies is not None else 0
n_live = 100 + 30 * len(lens_light_models) + 30 * n_extra + 30 * n_scaling
n_live = 100 + 30 * len(lens_dict) + 30 * n_extra + 30 * n_scaling

lens_dict = {f"lens_{i}": m for i, m in enumerate(lens_light_models)}
model = af.Collection(
galaxies=af.Collection(**lens_dict),
extra_galaxies=extra_galaxies,
Expand Down Expand Up @@ -266,24 +263,22 @@ def source_lp_1(

# --- main lens full models (light fixed from stage 0, mass + shear free) ---
# Only lens_0 carries the ExternalShear; one shear per group system.
lens_full_models = []
lens_dict = {}
for i in range(n_main):
lp0_lens = getattr(source_lp_result_0.instance.galaxies, f"lens_{i}")

mass = af.Model(al.mp.Isothermal)
mass.centre = lp0_lens.bulge.centre
mass.einstein_radius = af.UniformPrior(lower_limit=0.0, upper_limit=5.0)

lens_full_models.append(
af.Model(
al.Galaxy,
redshift=redshift_lens,
bulge=lp0_lens.bulge,
disk=lp0_lens.disk,
point=lp0_lens.point,
mass=mass,
shear=af.Model(al.mp.ExternalShear) if i == 0 else None,
)
lens_dict[f"lens_{i}"] = af.Model(
al.Galaxy,
redshift=redshift_lens,
bulge=lp0_lens.bulge,
disk=lp0_lens.disk,
point=lp0_lens.point,
mass=mass,
shear=af.Model(al.mp.ExternalShear) if i == 0 else None,
)

# --- extra lens galaxy models (light fixed, mass bounded by luminosity) ---
Expand Down Expand Up @@ -343,13 +338,12 @@ def source_lp_1(
af.Collection(scaling_mass_models) if scaling_mass_models else None
)

lens_dict = {f"lens_{i}": m for i, m in enumerate(lens_full_models)}
lens_dict["source"] = af.Model(
source = af.Model(
al.Galaxy, redshift=redshift_source, bulge=source_bulge
)

model = af.Collection(
galaxies=af.Collection(**lens_dict),
galaxies=af.Collection(**lens_dict, source=source),
extra_galaxies=extra_galaxies,
scaling_galaxies=scaling_galaxies,
)
Expand Down Expand Up @@ -473,7 +467,7 @@ def source_pix_1(
shear=lp_lens_model.shear,
)

lens_dict["source"] = af.Model(
source = af.Model(
al.Galaxy,
redshift=source_lp_result_1.instance.galaxies.source.redshift,
pixelization=af.Model(
Expand All @@ -486,7 +480,7 @@ def source_pix_1(
)

model = af.Collection(
galaxies=af.Collection(**lens_dict),
galaxies=af.Collection(**lens_dict, source=source),
extra_galaxies=source_lp_result_1.model.extra_galaxies,
scaling_galaxies=source_lp_result_1.model.scaling_galaxies,
)
Expand Down Expand Up @@ -598,7 +592,7 @@ def source_pix_2(
shear=pix1_lens_instance.shear,
)

lens_dict["source"] = af.Model(
source = af.Model(
al.Galaxy,
redshift=source_lp_result_1.instance.galaxies.source.redshift,
pixelization=af.Model(
Expand All @@ -611,7 +605,7 @@ def source_pix_2(
)

model = af.Collection(
galaxies=af.Collection(**lens_dict),
galaxies=af.Collection(**lens_dict, source=source),
extra_galaxies=source_pix_result_1.instance.extra_galaxies,
scaling_galaxies=source_pix_result_1.instance.scaling_galaxies,
)
Expand Down Expand Up @@ -710,10 +704,8 @@ def light_lp(
shear=lens_instance.shear,
)

lens_dict["source"] = source

model = af.Collection(
galaxies=af.Collection(**lens_dict),
galaxies=af.Collection(**lens_dict, source=source),
extra_galaxies=extra_galaxies,
scaling_galaxies=source_pix_result_2.instance.scaling_galaxies,
)
Expand Down Expand Up @@ -855,10 +847,8 @@ def mass_total(
shear=lens_model.shear,
)

lens_dict["source"] = source

model = af.Collection(
galaxies=af.Collection(**lens_dict),
galaxies=af.Collection(**lens_dict, source=source),
extra_galaxies=extra_galaxies,
scaling_galaxies=scaling_galaxies,
)
Expand Down
11 changes: 3 additions & 8 deletions scripts/group/start_here.py
Original file line number Diff line number Diff line change
Expand Up @@ -195,7 +195,7 @@
"""
# Main Lens Galaxies:

lens_models = []
lens_dict = {}

for i, centre in enumerate(main_lens_centres):

Expand All @@ -209,16 +209,14 @@
mass = af.Model(al.mp.Isothermal)
mass.centre = (centre[0], centre[1])

lens = af.Model(
lens_dict[f"lens_{i}"] = af.Model(
al.Galaxy,
redshift=0.5,
bulge=bulge,
mass=mass,
shear=af.Model(al.mp.ExternalShear) if i == 0 else None,
)

lens_models.append(lens)

# Source:

bulge = al.model_util.mge_model_from(
Expand All @@ -232,10 +230,7 @@

# Overall Lens Model:

lens_dict = {f"lens_{i}": m for i, m in enumerate(lens_models)}
lens_dict["source"] = source

model = af.Collection(galaxies=af.Collection(**lens_dict))
model = af.Collection(galaxies=af.Collection(**lens_dict, source=source))

"""
We can print the model to show the parameters that the model is composed of, which shows many of the MGE's fixed
Expand Down
45 changes: 34 additions & 11 deletions scripts/guides/advanced/add_a_profile.py
Original file line number Diff line number Diff line change
Expand Up @@ -111,14 +111,18 @@ def __init__(
self.slope = 2.0

@aa.grid_dec.to_vector_yx
@aa.grid_dec.transform
@aa.grid_dec.transform(rotate_back=True)
def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
"""
Calculate the deflection angles on a grid of (y,x) arc-second coordinates.

The input grid of (y,x) coordinates are transformed to a coordinate system centred on the profile centre with
and rotated based on the position angle defined from its `ell_comps` (this is described fully below).

Because this method computes deflection components using the rotated grid coordinates (i.e. the
components are expressed in the profile's frame), ``rotate_back=True`` is set so the decorator
automatically rotates them back to the observer frame.

The numerical backend can be selected via the ``xp`` argument, allowing this
method to be used with both NumPy and JAX (e.g. inside ``jax.jit``-compiled
code). This is described fully later in this example.
Expand Down Expand Up @@ -154,11 +158,7 @@ def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
psi,
)
)
return self.rotated_grid_from_reference_frame_from(
grid=xp.multiply(factor, xp.vstack((deflection_y, deflection_x)).T),
xp=xp,
**kwargs,
)
return xp.multiply(factor, xp.vstack((deflection_y, deflection_x)).T)


"""
Expand Down Expand Up @@ -446,11 +446,29 @@ def deflections_yx_2d_from(self, grid: aa.Grid2D, xp=np, **kwargs):
"""
__Do I Rotate Back?__

The ``aa.grid_dec.transform`` decorator does not rotate the deflection angles back to the original reference frame.
Whether the returned result needs back-rotation depends on what kind of quantity the function returns
and which coordinate frame its components are expressed in.

**Scalars** (convergence, potential): frame-invariant — no back-rotation needed.

**Vectors** (deflection angles): it depends on how the function computes them. The ``@aa.grid_dec.transform``
decorator transforms the input grid into the profile's rotated reference frame. If the function then computes
deflection components ``(alpha_y, alpha_x)`` using that rotated grid — as the ``Isothermal`` example above does —
those components are expressed in the profile's frame. Since the ray-tracing code expects observer-frame
components, they must be rotated back. Setting ``rotate_back=True`` handles this automatically by calling
``self.rotated_grid_from_reference_frame_from(...)`` on the result.

This is because deflection angles are typically applied in the same reference frame as the mass profile itself,
making such a rotation unnecessary. When implementing a custom mass profile, you should therefore ensure that any
required rotation of the deflection angles is correctly accounted for after they are calculated.
However, back-rotation is **not** always needed for vectors. If a function computes a scalar quantity in the
profile frame (e.g. a radial deflection magnitude) and then reconstructs the Cartesian vector using
observer-frame geometry, the result is already in the observer frame. In that case, ``rotate_back`` should
remain ``False``.

The rule is: look at which coordinate basis the returned components are expressed in. If they use the rotated
basis, set ``rotate_back=True``. If they are already observer-frame components, leave it as ``False``.

**Spin-2 quantities** (shear): these transform under a coordinate rotation by twice the profile angle (the
spin-2 transformation law). This is not handled by ``rotate_back`` — shear methods must apply the ``2 * angle``
rotation manually via ``self.rotated_grid_from_reference_frame_from(grid=..., angle=self.angle(xp) * 2)``.

__Lens Modeling__

Expand Down Expand Up @@ -579,10 +597,15 @@ def __init__(
# super().__init__(centre=centre, ell_comps=(0.0, 0.0))

@aa.grid_dec.to_vector_yx
@aa.grid_dec.transform
@aa.grid_dec.transform(rotate_back=True)
def deflections_yx_2d_from(self, grid: aa.type.Grid2DLike, xp=np, **kwargs):
"""
REQUIRED: The function is key for all lensing calculations and must be implemented.

Set ``rotate_back=True`` if your function computes deflection components using the rotated
grid coordinates (i.e. the components are expressed in the profile's frame). The decorator
will rotate them back to the observer frame automatically. If your function reconstructs
observer-frame components from scalar quantities, leave ``rotate_back=False``.
"""
pass

Expand Down
Loading