diff --git a/dataset/cluster/a2744/mass.csv b/dataset/cluster/a2744/mass.csv index b2d1e9c97..597b90f53 100644 --- a/dataset/cluster/a2744/mass.csv +++ b/dataset/cluster/a2744/mass.csv @@ -1,4 +1,4 @@ -galaxy,attr_name,profile_class,y,x,ra,rs,b0,mass_at_200,redshift_object,redshift_source,redshift -lens_0,mass,dPIEMassSph,-30.7476,14.0880,8.0,20.0,3.0,,,,0.308 -lens_1,mass,dPIEMassSph,-10.6992,-3.8597,8.0,20.0,3.0,,,,0.308 -host_halo,dark,NFWMCRLudlowSph,-30.7476,14.0880,,,,2000000000000000.0,0.308,5.662,0.308 +galaxy,attr_name,profile_class,y,x,sigma,r_core,r_cut,mass_at_200,redshift_object,redshift_source,H0,Om0,redshift +lens_0,mass,dPIEMassSph,-30.7476,14.0880,290.0,8.0,20.0,,0.308,5.662,67.66,0.30966,0.308 +lens_1,mass,dPIEMassSph,-10.6992,-3.8597,290.0,8.0,20.0,,0.308,5.662,67.66,0.30966,0.308 +host_halo,dark,NFWMCRLudlowSph,-30.7476,14.0880,,,,2000000000000000.0,0.308,5.662,,,0.308 diff --git a/dataset/cluster/a2744/prep.py b/dataset/cluster/a2744/prep.py index 6be237333..1d4ba8125 100644 --- a/dataset/cluster/a2744/prep.py +++ b/dataset/cluster/a2744/prep.py @@ -104,17 +104,22 @@ def fetch_tsv(table): with open(OUT / "mass.csv", "w", newline="") as f: writer = csv.writer(f) writer.writerow( - ["galaxy", "attr_name", "profile_class", "y", "x", "ra", "rs", "b0", - "mass_at_200", "redshift_object", "redshift_source", "redshift"] + ["galaxy", "attr_name", "profile_class", "y", "x", "sigma", "r_core", "r_cut", + "mass_at_200", "redshift_object", "redshift_source", "H0", "Om0", "redshift"] ) + # sigma is Lenstool's fiducial v_disp (km/s); 290 km/s corresponds to a lens + # strength b0 ~ 3.0" at (z_l = 0.308, z_s = 5.662). These are initialization + # placeholders — the modeling script promotes sigma/r_core/r_cut to priors. + # H0 / Om0 are the (fixed) cosmology constants entering the b0 normalization; + # they must be present so they load as fixed values rather than free priors. for i, bcg in enumerate((bcg_1, bcg_2)): writer.writerow( [f"lens_{i}", "mass", "dPIEMassSph", f"{bcg[1]:.4f}", f"{bcg[2]:.4f}", - 8.0, 20.0, 3.0, "", "", "", 0.308] + 290.0, 8.0, 20.0, "", 0.308, redshift_source_max, 67.66, 0.30966, 0.308] ) writer.writerow( ["host_halo", "dark", "NFWMCRLudlowSph", f"{bcg_1[1]:.4f}", f"{bcg_1[2]:.4f}", - "", "", "", 2e15, 0.308, redshift_source_max, 0.308] + "", "", "", 2e15, 0.308, redshift_source_max, "", "", 0.308] ) with open(OUT / "point.csv", "w", newline="") as f: diff --git a/dataset/cluster/simple/data.fits b/dataset/cluster/simple/data.fits index 134f04494..03049747a 100644 Binary files a/dataset/cluster/simple/data.fits and b/dataset/cluster/simple/data.fits differ diff --git a/dataset/cluster/simple/mass.csv b/dataset/cluster/simple/mass.csv index de830fba1..eefc820da 100644 --- a/dataset/cluster/simple/mass.csv +++ b/dataset/cluster/simple/mass.csv @@ -1,4 +1,4 @@ -galaxy,attr_name,profile_class,y,x,ra,rs,b0,mass_at_200,redshift_object,redshift_source,redshift -lens_0,mass,dPIEMassSph,0.0,0.0,8.0,20.0,3.0,,,,0.5 -lens_1,mass,dPIEMassSph,10.0,8.0,5.0,12.0,1.2,,,,0.5 -host_halo,dark,NFWMCRLudlowSph,0.0,0.0,,,,1995262314968882.8,0.5,2.0,0.5 +galaxy,attr_name,profile_class,y,x,sigma,r_core,r_cut,redshift_object,redshift_source,H0,Om0,mass_at_200,redshift +lens_0,mass,dPIEMassSph,0.0,0.0,330.0,8.0,20.0,0.5,2.0,67.66,0.30966,,0.5 +lens_1,mass,dPIEMassSph,10.0,8.0,210.0,5.0,12.0,0.5,2.0,67.66,0.30966,,0.5 +host_halo,dark,NFWMCRLudlowSph,0.0,0.0,,,,0.5,2.0,,,1995262314968882.8,0.5 diff --git a/dataset/cluster/simple/noise_map.fits b/dataset/cluster/simple/noise_map.fits index aa35403a7..210f06efb 100644 Binary files a/dataset/cluster/simple/noise_map.fits and b/dataset/cluster/simple/noise_map.fits differ diff --git a/dataset/cluster/simple/point_dataset_0.json b/dataset/cluster/simple/point_dataset_0.json index d7fc0cb44..8f60735a7 100644 --- a/dataset/cluster/simple/point_dataset_0.json +++ b/dataset/cluster/simple/point_dataset_0.json @@ -2,50 +2,50 @@ "type": "instance", "class_path": "autolens.point.dataset.PointDataset", "arguments": { - "time_delays_noise_map": null, - "redshift": 1.0, - "positions_noise_map": { + "time_delays": null, + "fluxes_noise_map": null, + "positions": { "type": "instance", - "class_path": "autoarray.structures.arrays.irregular.ArrayIrregular", + "class_path": "autoarray.structures.grids.irregular_2d.Grid2DIrregular", "arguments": { "values": { "type": "ndarray", "array": [ - 0.005, - 0.005, - 0.005 + [ + -9.15390625, + -19.291166924196176 + ], + [ + 0.03671875, + -0.07352194834211641 + ], + [ + 1.7234375000000002, + 23.35201416912914 + ] ], "dtype": "float64" } } }, - "fluxes_noise_map": null, + "name": "point_0", "fluxes": null, - "positions": { + "time_delays_noise_map": null, + "positions_noise_map": { "type": "instance", - "class_path": "autoarray.structures.grids.irregular_2d.Grid2DIrregular", + "class_path": "autoarray.structures.arrays.irregular.ArrayIrregular", "arguments": { "values": { "type": "ndarray", "array": [ - [ - -9.16875, - -19.276282112568627 - ], - [ - 0.0359375, - -0.0739730032399208 - ], - [ - 1.8125, - 23.341188851581837 - ] + 0.005, + 0.005, + 0.005 ], "dtype": "float64" } } }, - "time_delays": null, - "name": "point_0" + "redshift": 1.0 } } \ No newline at end of file diff --git a/dataset/cluster/simple/point_dataset_1.json b/dataset/cluster/simple/point_dataset_1.json index dfef8b699..b2661e898 100644 --- a/dataset/cluster/simple/point_dataset_1.json +++ b/dataset/cluster/simple/point_dataset_1.json @@ -2,50 +2,50 @@ "type": "instance", "class_path": "autolens.point.dataset.PointDataset", "arguments": { - "time_delays_noise_map": null, - "redshift": 2.0, - "positions_noise_map": { + "time_delays": null, + "fluxes_noise_map": null, + "positions": { "type": "instance", - "class_path": "autoarray.structures.arrays.irregular.ArrayIrregular", + "class_path": "autoarray.structures.grids.irregular_2d.Grid2DIrregular", "arguments": { "values": { "type": "ndarray", "array": [ - 0.005, - 0.005, - 0.005 + [ + -15.98203125, + 18.998883350418925 + ], + [ + 0.684375, + -0.5133004737014016 + ], + [ + 13.72578125, + -14.291674436932258 + ] ], "dtype": "float64" } } }, - "fluxes_noise_map": null, + "name": "point_1", "fluxes": null, - "positions": { + "time_delays_noise_map": null, + "positions_noise_map": { "type": "instance", - "class_path": "autoarray.structures.grids.irregular_2d.Grid2DIrregular", + "class_path": "autoarray.structures.arrays.irregular.ArrayIrregular", "arguments": { "values": { "type": "ndarray", "array": [ - [ - -15.96796875, - 19.00429600919258 - ], - [ - 0.68125, - -0.5133004737014016 - ], - [ - 13.69296875, - -14.308814523048826 - ] + 0.005, + 0.005, + 0.005 ], "dtype": "float64" } } }, - "time_delays": null, - "name": "point_1" + "redshift": 2.0 } } \ No newline at end of file diff --git a/dataset/cluster/simple/point_datasets.csv b/dataset/cluster/simple/point_datasets.csv index 10098f56f..c84814c8f 100644 --- a/dataset/cluster/simple/point_datasets.csv +++ b/dataset/cluster/simple/point_datasets.csv @@ -1,7 +1,7 @@ name,y,x,positions_noise,redshift -point_0,-9.16875,-19.276282112568627,0.005,1.0 -point_0,0.0359375,-0.0739730032399208,0.005,1.0 -point_0,1.8125,23.341188851581837,0.005,1.0 -point_1,-15.96796875,19.00429600919258,0.005,2.0 -point_1,0.68125,-0.5133004737014016,0.005,2.0 -point_1,13.69296875,-14.308814523048826,0.005,2.0 +point_0,-9.15390625,-19.291166924196176,0.005,1.0 +point_0,0.03671875,-0.07352194834211641,0.005,1.0 +point_0,1.7234375000000002,23.35201416912914,0.005,1.0 +point_1,-15.98203125,18.998883350418925,0.005,2.0 +point_1,0.684375,-0.5133004737014016,0.005,2.0 +point_1,13.72578125,-14.291674436932258,0.005,2.0 diff --git a/dataset/cluster/simple/tracer.json b/dataset/cluster/simple/tracer.json index c5fda962e..60f6ad830 100644 --- a/dataset/cluster/simple/tracer.json +++ b/dataset/cluster/simple/tracer.json @@ -2,6 +2,11 @@ "type": "instance", "class_path": "autolens.lens.tracer.Tracer", "arguments": { + "cosmology": { + "type": "instance", + "class_path": "autogalaxy.cosmology.model.Planck15", + "arguments": {} + }, "galaxies": { "type": "list", "values": [ @@ -9,13 +14,13 @@ "type": "instance", "class_path": "autogalaxy.galaxy.galaxy.Galaxy", "arguments": { - "label": null, "redshift": 0.5, + "label": null, "bulge": { "type": "instance", "class_path": "autogalaxy.profiles.light.standard.sersic.SersicSph", "arguments": { - "effective_radius": 3.0, + "intensity": 1.5, "centre": { "type": "tuple", "values": [ @@ -23,14 +28,18 @@ 0.0 ] }, - "sersic_index": 4.0, - "intensity": 1.5 + "effective_radius": 3.0, + "sersic_index": 4.0 } }, "mass": { "type": "instance", "class_path": "autogalaxy.profiles.mass.total.dual_pseudo_isothermal_mass.dPIEMassSph", "arguments": { + "r_cut": 20.0, + "sigma": 330.0, + "H0": 67.66, + "Om0": 0.30966, "centre": { "type": "tuple", "values": [ @@ -38,9 +47,9 @@ 0.0 ] }, - "b0": 3.0, - "ra": 8.0, - "rs": 20.0 + "redshift_object": 0.5, + "r_core": 8.0, + "redshift_source": 2.0 } } } @@ -49,13 +58,13 @@ "type": "instance", "class_path": "autogalaxy.galaxy.galaxy.Galaxy", "arguments": { - "label": null, "redshift": 0.5, + "label": null, "bulge": { "type": "instance", "class_path": "autogalaxy.profiles.light.standard.sersic.SersicSph", "arguments": { - "effective_radius": 1.5, + "intensity": 0.8, "centre": { "type": "tuple", "values": [ @@ -63,14 +72,18 @@ 8.0 ] }, - "sersic_index": 3.5, - "intensity": 0.8 + "effective_radius": 1.5, + "sersic_index": 3.5 } }, "mass": { "type": "instance", "class_path": "autogalaxy.profiles.mass.total.dual_pseudo_isothermal_mass.dPIEMassSph", "arguments": { + "r_cut": 12.0, + "sigma": 210.0, + "H0": 67.66, + "Om0": 0.30966, "centre": { "type": "tuple", "values": [ @@ -78,9 +91,9 @@ 8.0 ] }, - "b0": 1.2, - "ra": 5.0, - "rs": 12.0 + "redshift_object": 0.5, + "r_core": 5.0, + "redshift_source": 2.0 } } } @@ -89,13 +102,13 @@ "type": "instance", "class_path": "autogalaxy.galaxy.galaxy.Galaxy", "arguments": { - "label": null, "redshift": 0.5, + "label": null, "bulge": { "type": "instance", "class_path": "autogalaxy.profiles.light.standard.sersic.SersicSph", "arguments": { - "effective_radius": 0.8, + "intensity": 0.4, "centre": { "type": "tuple", "values": [ @@ -103,14 +116,18 @@ -6.5 ] }, - "sersic_index": 3.0, - "intensity": 0.4 + "effective_radius": 0.8, + "sersic_index": 3.0 } }, "mass": { "type": "instance", "class_path": "autogalaxy.profiles.mass.total.dual_pseudo_isothermal_mass.dPIEMassSph", "arguments": { + "r_cut": 9.992797406132079, + "sigma": 67.5980119451993, + "H0": 67.66, + "Om0": 0.30966, "centre": { "type": "tuple", "values": [ @@ -118,9 +135,9 @@ -6.5 ] }, - "b0": 0.12016655108639841, - "ra": 0.09992797406132078, - "rs": 9.992797406132079 + "redshift_object": 0.5, + "r_core": 0.09992797406132078, + "redshift_source": 2.0 } } } @@ -129,13 +146,13 @@ "type": "instance", "class_path": "autogalaxy.galaxy.galaxy.Galaxy", "arguments": { - "label": null, "redshift": 0.5, + "label": null, "bulge": { "type": "instance", "class_path": "autogalaxy.profiles.light.standard.sersic.SersicSph", "arguments": { - "effective_radius": 0.8, + "intensity": 0.32, "centre": { "type": "tuple", "values": [ @@ -143,14 +160,18 @@ 3.0 ] }, - "sersic_index": 3.0, - "intensity": 0.32 + "effective_radius": 0.8, + "sersic_index": 3.0 } }, "mass": { "type": "instance", "class_path": "autogalaxy.profiles.mass.total.dual_pseudo_isothermal_mass.dPIEMassSph", "arguments": { + "r_cut": 8.937829714197962, + "sigma": 63.9302525824687, + "H0": 67.66, + "Om0": 0.30966, "centre": { "type": "tuple", "values": [ @@ -158,9 +179,9 @@ 3.0 ] }, - "b0": 0.10748023074035522, - "ra": 0.0893782971419796, - "rs": 8.937829714197962 + "redshift_object": 0.5, + "r_core": 0.0893782971419796, + "redshift_source": 2.0 } } } @@ -169,13 +190,13 @@ "type": "instance", "class_path": "autogalaxy.galaxy.galaxy.Galaxy", "arguments": { - "label": null, "redshift": 0.5, + "label": null, "bulge": { "type": "instance", "class_path": "autogalaxy.profiles.light.standard.sersic.SersicSph", "arguments": { - "effective_radius": 0.8, + "intensity": 0.25, "centre": { "type": "tuple", "values": [ @@ -183,14 +204,18 @@ -5.0 ] }, - "sersic_index": 3.0, - "intensity": 0.25 + "effective_radius": 0.8, + "sersic_index": 3.0 } }, "mass": { "type": "instance", "class_path": "autogalaxy.profiles.mass.total.dual_pseudo_isothermal_mass.dPIEMassSph", "arguments": { + "r_cut": 7.9, + "sigma": 60.10407640085654, + "H0": 67.66, + "Om0": 0.30966, "centre": { "type": "tuple", "values": [ @@ -198,9 +223,9 @@ -5.0 ] }, - "b0": 0.095, - "ra": 0.079, - "rs": 7.9 + "redshift_object": 0.5, + "r_core": 0.079, + "redshift_source": 2.0 } } } @@ -209,13 +234,13 @@ "type": "instance", "class_path": "autogalaxy.galaxy.galaxy.Galaxy", "arguments": { - "label": null, "redshift": 0.5, + "label": null, "bulge": { "type": "instance", "class_path": "autogalaxy.profiles.light.standard.sersic.SersicSph", "arguments": { - "effective_radius": 0.8, + "intensity": 0.2, "centre": { "type": "tuple", "values": [ @@ -223,14 +248,18 @@ -9.0 ] }, - "sersic_index": 3.0, - "intensity": 0.2 + "effective_radius": 0.8, + "sersic_index": 3.0 } }, "mass": { "type": "instance", "class_path": "autogalaxy.profiles.mass.total.dual_pseudo_isothermal_mass.dPIEMassSph", "arguments": { + "r_cut": 7.065974808899336, + "sigma": 56.84292592299587, + "H0": 67.66, + "Om0": 0.30966, "centre": { "type": "tuple", "values": [ @@ -238,9 +267,9 @@ -9.0 ] }, - "b0": 0.084970583144992, - "ra": 0.07065974808899335, - "rs": 7.065974808899336 + "redshift_object": 0.5, + "r_core": 0.07065974808899335, + "redshift_source": 2.0 } } } @@ -249,13 +278,13 @@ "type": "instance", "class_path": "autogalaxy.galaxy.galaxy.Galaxy", "arguments": { - "label": null, "redshift": 0.5, + "label": null, "bulge": { "type": "instance", "class_path": "autogalaxy.profiles.light.standard.sersic.SersicSph", "arguments": { - "effective_radius": 0.8, + "intensity": 0.16, "centre": { "type": "tuple", "values": [ @@ -263,14 +292,18 @@ 13.0 ] }, - "sersic_index": 3.0, - "intensity": 0.16 + "effective_radius": 0.8, + "sersic_index": 3.0 } }, "mass": { "type": "instance", "class_path": "autogalaxy.profiles.mass.total.dual_pseudo_isothermal_mass.dPIEMassSph", "arguments": { + "r_cut": 6.32, + "sigma": 53.75872022286245, + "H0": 67.66, + "Om0": 0.30966, "centre": { "type": "tuple", "values": [ @@ -278,9 +311,9 @@ 13.0 ] }, - "b0": 0.07600000000000001, - "ra": 0.0632, - "rs": 6.32 + "redshift_object": 0.5, + "r_core": 0.0632, + "redshift_source": 2.0 } } } @@ -289,13 +322,13 @@ "type": "instance", "class_path": "autogalaxy.galaxy.galaxy.Galaxy", "arguments": { - "label": null, "redshift": 0.5, + "label": null, "bulge": { "type": "instance", "class_path": "autogalaxy.profiles.light.standard.sersic.SersicSph", "arguments": { - "effective_radius": 0.8, + "intensity": 0.13, "centre": { "type": "tuple", "values": [ @@ -303,14 +336,18 @@ 4.0 ] }, - "sersic_index": 3.0, - "intensity": 0.13 + "effective_radius": 0.8, + "sersic_index": 3.0 } }, "mass": { "type": "instance", "class_path": "autogalaxy.profiles.mass.total.dual_pseudo_isothermal_mass.dPIEMassSph", "arguments": { + "r_cut": 5.696771015233104, + "sigma": 51.03930638755519, + "H0": 67.66, + "Om0": 0.30966, "centre": { "type": "tuple", "values": [ @@ -318,9 +355,9 @@ 4.0 ] }, - "b0": 0.0685054742338158, - "ra": 0.056967710152331034, - "rs": 5.696771015233104 + "redshift_object": 0.5, + "r_core": 0.056967710152331034, + "redshift_source": 2.0 } } } @@ -329,13 +366,13 @@ "type": "instance", "class_path": "autogalaxy.galaxy.galaxy.Galaxy", "arguments": { - "label": null, "redshift": 0.5, + "label": null, "bulge": { "type": "instance", "class_path": "autogalaxy.profiles.light.standard.sersic.SersicSph", "arguments": { - "effective_radius": 0.8, + "intensity": 0.1, "centre": { "type": "tuple", "values": [ @@ -343,14 +380,18 @@ 9.0 ] }, - "sersic_index": 3.0, - "intensity": 0.1 + "effective_radius": 0.8, + "sersic_index": 3.0 } }, "mass": { "type": "instance", "class_path": "autogalaxy.profiles.mass.total.dual_pseudo_isothermal_mass.dPIEMassSph", "arguments": { + "r_cut": 4.996398703066039, + "sigma": 47.79901264117967, + "H0": 67.66, + "Om0": 0.30966, "centre": { "type": "tuple", "values": [ @@ -358,9 +399,9 @@ 9.0 ] }, - "b0": 0.060083275543199206, - "ra": 0.04996398703066039, - "rs": 4.996398703066039 + "redshift_object": 0.5, + "r_core": 0.04996398703066039, + "redshift_source": 2.0 } } } @@ -369,13 +410,13 @@ "type": "instance", "class_path": "autogalaxy.galaxy.galaxy.Galaxy", "arguments": { - "label": null, "redshift": 0.5, + "label": null, "bulge": { "type": "instance", "class_path": "autogalaxy.profiles.light.standard.sersic.SersicSph", "arguments": { - "effective_radius": 0.8, + "intensity": 0.08, "centre": { "type": "tuple", "values": [ @@ -383,14 +424,18 @@ -12.0 ] }, - "sersic_index": 3.0, - "intensity": 0.08 + "effective_radius": 0.8, + "sersic_index": 3.0 } }, "mass": { "type": "instance", "class_path": "autogalaxy.profiles.mass.total.dual_pseudo_isothermal_mass.dPIEMassSph", "arguments": { + "r_cut": 4.468914857098981, + "sigma": 45.20551512403241, + "H0": 67.66, + "Om0": 0.30966, "centre": { "type": "tuple", "values": [ @@ -398,9 +443,9 @@ -12.0 ] }, - "b0": 0.05374011537017761, - "ra": 0.0446891485709898, - "rs": 4.468914857098981 + "redshift_object": 0.5, + "r_core": 0.0446891485709898, + "redshift_source": 2.0 } } } @@ -409,13 +454,13 @@ "type": "instance", "class_path": "autogalaxy.galaxy.galaxy.Galaxy", "arguments": { - "label": null, "redshift": 0.5, + "label": null, "bulge": { "type": "instance", "class_path": "autogalaxy.profiles.light.standard.sersic.SersicSph", "arguments": { - "effective_radius": 0.8, + "intensity": 0.06, "centre": { "type": "tuple", "values": [ @@ -423,14 +468,18 @@ 5.5 ] }, - "sersic_index": 3.0, - "intensity": 0.06 + "effective_radius": 0.8, + "sersic_index": 3.0 } }, "mass": { "type": "instance", "class_path": "autogalaxy.profiles.mass.total.dual_pseudo_isothermal_mass.dPIEMassSph", "arguments": { + "r_cut": 3.8701937935974215, + "sigma": 42.068472032638006, + "H0": 67.66, + "Om0": 0.30966, "centre": { "type": "tuple", "values": [ @@ -438,9 +487,9 @@ 5.5 ] }, - "b0": 0.04654030511288038, - "ra": 0.038701937935974214, - "rs": 3.8701937935974215 + "redshift_object": 0.5, + "r_core": 0.038701937935974214, + "redshift_source": 2.0 } } } @@ -449,13 +498,13 @@ "type": "instance", "class_path": "autogalaxy.galaxy.galaxy.Galaxy", "arguments": { - "label": null, "redshift": 0.5, + "label": null, "bulge": { "type": "instance", "class_path": "autogalaxy.profiles.light.standard.sersic.SersicSph", "arguments": { - "effective_radius": 0.8, + "intensity": 0.05, "centre": { "type": "tuple", "values": [ @@ -463,14 +512,18 @@ 11.0 ] }, - "sersic_index": 3.0, - "intensity": 0.05 + "effective_radius": 0.8, + "sersic_index": 3.0 } }, "mass": { "type": "instance", "class_path": "autogalaxy.profiles.mass.total.dual_pseudo_isothermal_mass.dPIEMassSph", "arguments": { + "r_cut": 3.532987404449668, + "sigma": 40.19401838263497, + "H0": 67.66, + "Om0": 0.30966, "centre": { "type": "tuple", "values": [ @@ -478,9 +531,9 @@ 11.0 ] }, - "b0": 0.042485291572496, - "ra": 0.035329874044496674, - "rs": 3.532987404449668 + "redshift_object": 0.5, + "r_core": 0.035329874044496674, + "redshift_source": 2.0 } } } @@ -489,12 +542,14 @@ "type": "instance", "class_path": "autogalaxy.galaxy.galaxy.Galaxy", "arguments": { - "label": null, "redshift": 0.5, + "label": null, "dark": { "type": "instance", "class_path": "autogalaxy.profiles.mass.dark.nfw_mcr.NFWMCRLudlowSph", "arguments": { + "redshift_source": 2.0, + "mass_at_200": 1995262314968882.8, "centre": { "type": "tuple", "values": [ @@ -502,8 +557,6 @@ 0.0 ] }, - "mass_at_200": 1995262314968882.8, - "redshift_source": 2.0, "redshift_object": 0.5 } } @@ -513,20 +566,13 @@ "type": "instance", "class_path": "autogalaxy.galaxy.galaxy.Galaxy", "arguments": { - "label": null, "redshift": 1.0, + "label": null, "bulge": { "type": "instance", "class_path": "autogalaxy.profiles.light.standard.sersic_core.SersicCore", "arguments": { - "effective_radius": 0.3, - "centre": { - "type": "tuple", - "values": [ - 0.3, - 0.5 - ] - }, + "intensity": 2.0, "ell_comps": { "type": "tuple", "values": [ @@ -535,10 +581,17 @@ ] }, "alpha": 3.0, + "centre": { + "type": "tuple", + "values": [ + 0.3, + 0.5 + ] + }, + "gamma": 0.25, "radius_break": 0.025, "sersic_index": 1.0, - "intensity": 2.0, - "gamma": 0.25 + "effective_radius": 0.3 } } } @@ -547,20 +600,13 @@ "type": "instance", "class_path": "autogalaxy.galaxy.galaxy.Galaxy", "arguments": { - "label": null, "redshift": 2.0, + "label": null, "bulge": { "type": "instance", "class_path": "autogalaxy.profiles.light.standard.sersic_core.SersicCore", "arguments": { - "effective_radius": 0.3, - "centre": { - "type": "tuple", - "values": [ - -0.8, - 1.2 - ] - }, + "intensity": 2.0, "ell_comps": { "type": "tuple", "values": [ @@ -569,20 +615,22 @@ ] }, "alpha": 3.0, + "centre": { + "type": "tuple", + "values": [ + -0.8, + 1.2 + ] + }, + "gamma": 0.25, "radius_break": 0.025, "sersic_index": 1.0, - "intensity": 2.0, - "gamma": 0.25 + "effective_radius": 0.3 } } } } ] - }, - "cosmology": { - "type": "instance", - "class_path": "autogalaxy.cosmology.model.Planck15", - "arguments": {} } } } \ No newline at end of file diff --git a/llms-full.txt b/llms-full.txt index c1d0ac4fe..737afdbdc 100644 --- a/llms-full.txt +++ b/llms-full.txt @@ -137,6 +137,8 @@ AUTO-GENERATED by PyAutoBuild — do not edit by hand; regenerate with generate. - Contents: Prerequisites, SOURCE LP PIPELINE, SOURCE PIX PIPELINE 1, SOURCE PIX PIPELINE 2, LIGHT LP PIPELINE, MASS TOTAL PIPELINE, Dataset, Settings AutoFit, Redshifts, Mesh Shape, SLaM Pipeline - [Pixelization: Source Science](scripts/imaging/features/pixelization/source_science.py): Source science focuses on studying the highly magnified properties of the background lensed source galaxy (or galaxies). - Contents: Model Fit, Interpolated Source, Source Flux, Zoom, Errors, Magnification, Masking, Magnification via Mesh, Reconstruction CSV +- [__Log Likelihood Function: Potential Correction (Gravitational Imaging)__](scripts/imaging/features/potential_correction/likelihood_function.py): This script provides a step-by-step guide of the **PyAutoLens** potential-correction `log_likelihood_function` (the Bayesian evidence of `al.pc.FitDpsiSrcImaging`), which jointly inverts `Imaging` data for a pixelized source and pixelized corrections $\delta\psi$ to the lensing potential. + - Contents: Simplifications, Prerequisites, Dataset, Arc Mask, Smooth Starting Model, Dpsi Mesh, Dpsi Gradient Operator, Source Gradients, PSF Blur Matrix, Dpsi Mapping Matrix, Source Inversion Blocks, Joint System, Solve, Dkappa Map, Evidence Terms, Verification, JAX / xp, Wrap Up - [Features: Scaling Relation Fit](scripts/imaging/features/scaling_relation/fit.py): A strong lens system often has many foreground galaxies near the line of sight to the source, in addition to the primary lens. As the number of foreground galaxies grows, modelling each one individually with its own free `einstein_radius` parameter rapidly becomes intractable. - Contents: Prerequisites, Dataset & Mask, Over Sampling, Centres + Luminosities, MGE Basis, Galaxies, Tracer, Scaling Relation Tour, Intensities, Wrap Up - [__Log Likelihood Function: Scaling Relation__](scripts/imaging/features/scaling_relation/likelihood_function.py): This script describes the additional steps required to compute the `log_likelihood` for a strong lens whose foreground galaxy population is split between two tiers — individually-modelled extras (each with its own free `einstein_radius`) and scaling-tier extras (whose Einstein radii are derived from a shared reference-anchored relation `einstein_radius = einstein_radius_ref * (luminosity / reference_luminosity) ** 0.5`). @@ -331,6 +333,9 @@ AUTO-GENERATED by PyAutoBuild — do not edit by hand; regenerate with generate. - Contents: Prerequisites, Interferometer SLaM Description, High Resolution Dataset, SOURCE LP PIPELINE, SOURCE PIX PIPELINE 1, SOURCE PIX PIPELINE 2, MASS TOTAL PIPELINE, Two Datasets, Sparse Operators, Settings, Settings AutoFit, Redshifts, Mesh Shape, SLaM Pipeline - [Pixelization: Source Reconstruction](scripts/interferometer/features/pixelization/source_science.py): Source science focuses on studying the highly magnified properties of the background lensed source galaxy (or galaxies). - Contents: Model Fit, Interpolated Source, Source Flux, Zoom, Errors, Magnification, Masking, Magnification via Mesh, Reconstruction CSV +- [__Log Likelihood Function: Potential Correction (Interferometer)__](scripts/interferometer/features/potential_correction/likelihood_function.py): This script provides a step-by-step guide of the **PyAutoLens** visibility-space potential-correction `log_likelihood_function` (the Bayesian evidence of `al.pc.FitDpsiSrcInterferometer`), which jointly inverts `Interferometer` data for a pixelized source and pixelized corrections $\delta\psi$ to the lensing potential. + - Contents: Simplifications, Prerequisites, Dataset, Source Inversion Blocks, Dpsi Mesh, Correction Response, Joint Response, Dense Route, Sparse (w-tilde) Route, Solve, Dkappa Map, Evidence Terms, Verification, JAX / xp, Wrap Up +- [Feature: Potential Correction (Gravitational Imaging)](scripts/interferometer/features/potential_correction/start_here.py): This example performs gravitational imaging on interferometer data: pixelized corrections $\delta\psi$ to the lensing potential are reconstructed jointly with the pixelized source, revealing mass structure (e.g. dark subhaloes) the smooth lens model omits, as a convergence-correction map $\delta\kappa = \frac{1}{2}\nabla^2\delta\psi$. - [Subhalo Detection: Start Here](scripts/interferometer/features/subhalo/detect/start_here.py): Strong gravitational lenses can be used to detect the presence of small-scale dark matter (DM) subhalos. This occurs when the DM subhalo overlaps the lensed source emission, and therefore gravitationally perturbs the observed image of the lensed source galaxy. - Contents: SLaM Pipelines, Grid Search, Pixelized Source, Model, SOURCE LP PIPELINE, SOURCE PIX PIPELINE 1, SOURCE PIX PIPELINE 2, MASS TOTAL PIPELINE, Two Datasets, Sparse Operators, Settings, Settings AutoFit, Redshifts, Mesh Shape, SLaM Pipeline - [start_here](scripts/interferometer/features/subhalo/sensitivity/start_here.py): (no summary in script docstring) @@ -371,7 +376,7 @@ AUTO-GENERATED by PyAutoBuild — do not edit by hand; regenerate with generate. - Contents: Model, Dataset Paths, Ray Tracing, Point Solver, Point Datasets, Visualize, Tracer json, Imaging, Fluxes, Point Dataset, Time Delays - [Simulator: Hubble Constant Time Delays](scripts/point_source/simulator_sample.py): A multiply imaged lensed quasar with time delays can measure the Hubble constant, which is a fundamental Cosmological parameter that describes the rate of expansion of the universe. This is because the the difference between the geometric time delay and the physical time delay is proportional to the Hubble constant. - Contents: Model, Dataset Paths, Point Solver, Sample Model Distributions, Simulate, Sample Instances -- [Start Here: Imaging](scripts/point_source/start_here.py): Strong gravitational lenses often have point sources (e.g. quasars) that are being lensed, appearing as two or four distinct point-like images. These lenses are particularly useful for measuring cosmological parameters like the Hubble constant, and for studying the small-scale properties of dark matter. +- [Start Here: Point Source](scripts/point_source/start_here.py): Strong gravitational lenses often have point sources (e.g. quasars) that are being lensed, appearing as two or four distinct point-like images. These lenses are particularly useful for measuring cosmological parameters like the Hubble constant, and for studying the small-scale properties of dark matter. - Contents: JAX, Google Colab Setup, Imports, Dataset, Point Solver, Model, Name Pairing, Model Fit, Live Visual Update, Result, Model Your Own Lens, Fluxes and Time Delays, Simulator, Sample, Wrap Up ## cluster @@ -401,6 +406,8 @@ AUTO-GENERATED by PyAutoBuild — do not edit by hand; regenerate with generate. - Contents: Benefits, Interferometer Masking, Interferometer Dataset, Imaging Dataset, Imaging Masking, Analysis, Model, Shared Source Mesh (Pixelization), Search, Result - [Simulator: Mutli Interferometer](scripts/multi/features/imaging_and_interferometer/simulator.py): This script simulates `Interferometer` data of a 'galaxy-scale' strong lens where: - Contents: Simulate, Ray Tracing, Output, Visualize, Tracer json +- [Modeling: Combined Imaging + Point Source](scripts/multi/features/imaging_and_point_source/modeling.py): This script jointly fits CCD imaging of a lensed quasar's **extended arcs** and its **point-source observables** (image positions and time delays) with a **single lens mass model**, using PyAutoFit's factor-graph API to sample the joint likelihood with one non-linear search. + - Contents: Dataset, Quasar Image Masking, Model, Analysis Factor & Factor Graph, Search & Model-Fit, Result - [Modeling Features: One By One](scripts/multi/features/one_by_one/modeling.py): Multi-wavelength analysis does not necessarily require us to fit all datasets simultaneously. Instead, we can fit one dataset first in order to infer a robust lens and source model, and then fit the next dataset, using the inferred model as the starting point. - Contents: Model, Colors, Pixel Scales, Dataset & Mask, Analysis, Search, Result, Second Dataset Mass Model Fixed, Second Dataset Offset - [Modeling: Pixelized](scripts/multi/features/pixelization/modeling.py): This script fits a multi-wavelength `Imaging` dataset of a 'galaxy-scale' strong lens with a model where: @@ -430,6 +437,8 @@ AUTO-GENERATED by PyAutoBuild — do not edit by hand; regenerate with generate. ## weak +- [Real Data: Combined Strong + Weak Lensing of Abell 2744](scripts/weak/features/strong_lensing/a2744.py): This script constrains a **single cluster-scale mass model** of **Abell 2744** ("Pandora's Cluster", z = 0.308) using both its **strong-lensing** multiple images and its **weak-lensing** shear catalogue — real data on both sides, fitted jointly with one non-linear search via PyAutoFit's factor-graph API. + - Contents: Dataset (Strong), Dataset (Weak), Model, Analysis Factors & Factor Graph, Search & Model-Fit, Result - [Fit: Combined Strong + Weak Lensing](scripts/weak/features/strong_lensing/fit.py): This script fits the combined strong+weak dataset simulated by `simulator.py` in this folder with a single shared `Tracer`: the imaging data via `FitImaging` and the shear catalogue via `FitWeak`. - Contents: Dataset, Tracer, Imaging Fit, Weak Fit, Joint Likelihood, Shear Profile - [Modeling: Combined Strong + Weak Lensing](scripts/weak/features/strong_lensing/modeling.py): This script fits the combined dataset of `simulator.py` — an `Imaging` dataset of strongly lensed arcs and a `WeakDataset` of the surrounding shear field — with a **single lens mass model**, using PyAutoFit's factor-graph API to sample the joint likelihood with one non-linear search. @@ -446,6 +455,8 @@ AUTO-GENERATED by PyAutoBuild — do not edit by hand; regenerate with generate. - Contents: Catalogue Download, Catalogue Load & Projection, Quality Cuts, Weak Dataset, Mass Map, Model & Search, Result - [Simulator: Weak Lensing](scripts/weak/simulator.py): This script simulates a weak gravitational lensing shear catalogue. Unlike the imaging simulator (which produces a 2D image of the lensed source) the weak-lensing simulator produces a *catalogue* of (gamma_2, gamma_1) shear measurements at the (y, x) positions of a population of background source galaxies. - Contents: Model, Dataset Paths, Ray Tracing, Source Positions, Simulator, Output, Visualize +- [Start Here: Weak Lensing](scripts/weak/start_here.py): Gravitational lensing does not stop at the spectacular arcs and multiple images of strong lensing. Every galaxy behind a massive structure is slightly sheared by its gravity — a percent-level distortion invisible in any single galaxy but measurable statistically across thousands. This is **weak lensing**, and it traces mass on scales far beyond the strong-lensing core. + - Contents: JAX, Google Colab Setup, Imports, Catalogue Download, Catalogue Load & Projection, Quality Cuts, Weak Dataset, Mass Map, Model, Model Fit, Result, Model Your Own Catalogue, Simulator, Wrap Up ## guides @@ -459,6 +470,8 @@ AUTO-GENERATED by PyAutoBuild — do not edit by hand; regenerate with generate. - Contents: Illustration, Numerics, Images, Adaptive Over Sampling, Multiple Lens Galaxies, Ray Tracing, Default Ray Tracing, Dataset & Modeling, Pixelization - [Chaining: Over Sample](scripts/guides/advanced/over_sampling_chaining.py): Over sampling is a numerical technique where the images of light profiles and galaxies are evaluated on a higher resolution grid than the image data to ensure the calculation is accurate. - Contents: Start Here Notebook, Dataset + Masking, Paths, Redshifts, Model (Search 1), Search + Analysis + Model-Fit (Search 1), Result (Search 1), Over Sampling (Search 2) +- [Potential Corrections (Gravitational Imaging)](scripts/guides/advanced/potential_correction.py): Parametric lens mass models (e.g. the `PowerLaw` used throughout the workspace) are smooth: their convergence and deflection angles vary gradually across the image plane. Real lens galaxies are not perfectly smooth. They contain substructure — dark matter subhaloes, globular clusters, satellite galaxies — and their large-scale mass distributions can depart from ellipsoidal symmetry in ways no simple parametric profile captures. + - Contents: Simulate, Smooth Model Fit, Joint Fit, Dkappa Map, Iterative Fit, Evidence Sampling - [COOLEST](scripts/guides/coolest_interop.py): COOLEST (COde-independent Organized LEns STandard, https://github.com/aymgal/COOLEST) is a standard for storing and exchanging strong lens models between different lens modeling software, for example lenstronomy, herculens and GLEE. - Contents: Conventions, Lens Model, Export, Import, Round Trip, NFW Profiles - [Data Structures](scripts/guides/data_structures.py): This tutorial illustrates the data structure objects which data and results quantities are stored using, which are extensions of NumPy arrays. diff --git a/notebooks/cluster/csv_api.ipynb b/notebooks/cluster/csv_api.ipynb index 19638d6f0..d90e22626 100644 --- a/notebooks/cluster/csv_api.ipynb +++ b/notebooks/cluster/csv_api.ipynb @@ -293,8 +293,26 @@ "# family carries the 2 source-galaxy Point components.\n", "\n", "mass_profiles = {\n", - " \"lens_0\": {\"mass\": al.mp.dPIEMassSph(centre=(0.0, 0.0), ra=8.0, rs=20.0, b0=3.0)},\n", - " \"lens_1\": {\"mass\": al.mp.dPIEMassSph(centre=(10.0, 8.0), ra=5.0, rs=12.0, b0=1.2)},\n", + " \"lens_0\": {\n", + " \"mass\": al.mp.dPIEMassSph(\n", + " centre=(0.0, 0.0),\n", + " sigma=330.0,\n", + " r_core=8.0,\n", + " r_cut=20.0,\n", + " redshift_object=redshift_lens,\n", + " redshift_source=max(source_redshifts),\n", + " )\n", + " },\n", + " \"lens_1\": {\n", + " \"mass\": al.mp.dPIEMassSph(\n", + " centre=(10.0, 8.0),\n", + " sigma=210.0,\n", + " r_core=5.0,\n", + " r_cut=12.0,\n", + " redshift_object=redshift_lens,\n", + " redshift_source=max(source_redshifts),\n", + " )\n", + " },\n", " \"host_halo\": {\n", " \"dark\": al.mp.NFWMCRLudlowSph(\n", " centre=(0.0, 0.0),\n", @@ -485,9 +503,9 @@ "mutate the returned model:\n", "\n", "```python\n", - "galaxy_models[\"lens_0\"].mass.ra = af.UniformPrior(lower_limit=1.0, upper_limit=15.0)\n", - "galaxy_models[\"lens_0\"].mass.rs = af.UniformPrior(lower_limit=5.0, upper_limit=40.0)\n", - "galaxy_models[\"lens_0\"].mass.b0 = af.UniformPrior(lower_limit=0.1, upper_limit=10.0)\n", + "galaxy_models[\"lens_0\"].mass.sigma = af.UniformPrior(lower_limit=50.0, upper_limit=600.0)\n", + "galaxy_models[\"lens_0\"].mass.r_core = af.UniformPrior(lower_limit=1.0, upper_limit=15.0)\n", + "galaxy_models[\"lens_0\"].mass.r_cut = af.UniformPrior(lower_limit=5.0, upper_limit=40.0)\n", "```\n", "\n", "This is the same composition pattern as ``af.Model(al.Galaxy, mass=...)``\n", @@ -505,11 +523,12 @@ " point_table,\n", ")\n", "\n", - "# Mutate selected params on the main-lens mass profiles into priors.\n", + "# Mutate selected params on the main-lens mass profiles into priors (sigma is\n", + "# Lenstool's fiducial v_disp in km/s; radii in arcsec).\n", "for name in (\"lens_0\", \"lens_1\"):\n", - " galaxy_models[name].mass.ra = af.UniformPrior(lower_limit=1.0, upper_limit=15.0)\n", - " galaxy_models[name].mass.rs = af.UniformPrior(lower_limit=5.0, upper_limit=40.0)\n", - " galaxy_models[name].mass.b0 = af.UniformPrior(lower_limit=0.1, upper_limit=10.0)\n", + " galaxy_models[name].mass.sigma = af.UniformPrior(lower_limit=50.0, upper_limit=600.0)\n", + " galaxy_models[name].mass.r_core = af.UniformPrior(lower_limit=1.0, upper_limit=15.0)\n", + " galaxy_models[name].mass.r_cut = af.UniformPrior(lower_limit=5.0, upper_limit=40.0)\n", "\n", "# Host halo: free mass_at_200, keep centre + redshifts fixed.\n", "galaxy_models[\"host_halo\"].dark.mass_at_200 = af.LogUniformPrior(\n", @@ -538,7 +557,7 @@ "member and emitting an ``attr_name`` column would be more overhead than\n", "signal \u2014 every row uses the same ``dPIEMassSph`` mass profile with\n", "parameters derived from the reference-anchored scaling relation's shared\n", - "``b0_ref`` normalization (see ``modeling.py``).\n", + "``sigma_ref`` normalization (see ``modeling.py``).\n", "\n", "``al.galaxy_table_to_csv`` and ``al.galaxy_table_from_csv`` are the\n", "schema-specific writers/readers. The simulator emits 10 scaling members\n", @@ -596,18 +615,18 @@ "\n", "__Lenstool-Parameterized Rows__\n", "\n", - "Because the ``profile_class`` column dispatches against the full ``al.mp`` namespace, a\n", - "``mass.csv`` can carry rows in **Lenstool's native parameterization** via ``dPIEMassLenstool``\n", - "\u2014 the columns become the ``.par``-file keywords verbatim::\n", + "The default dPIE classes ARE Lenstool's native parameterization, so an elliptical ``dPIEMass``\n", + "row's columns are the ``.par``-file keywords verbatim::\n", "\n", " galaxy,attr_name,profile_class,y,x,ellipticity,angle_pos,sigma,r_core,r_cut,redshift_object,redshift_source,H0,Om0,redshift\n", - " O1,mass,dPIEMassLenstool,1.479,-2.997,0.678,8.971,987.34,18.96,283.54,0.39,11.76,70.0,0.3,0.39\n", + " O1,mass,dPIEMass,1.479,-2.997,0.678,8.971,987.34,18.96,283.54,0.39,11.76,70.0,0.3,0.39\n", "\n", "``sigma`` is Lenstool's fiducial ``v_disp`` (sigma_LT), radii are in arcsec, and the run's own\n", "cosmology travels as the flat ``H0`` / ``Om0`` columns. ``scripts/cluster/lenstool/`` builds its\n", "entire 149-component published model this way \u2014 the ``.par`` file becomes one canonical CSV. Note\n", "the multi-plane convention: ``redshift_source`` must be the tracer's *final* (highest) source\n", - "plane.\n", + "plane. (The internal ``(ra, rs, b0)`` parameterization remains available for CSV rows via the\n", + "non-standard ``dPIEMassB0`` / ``dPIEMassB0Sph`` classes.)\n", "\n", "Light-profile CSVs (``light.csv``) support the linear / operated variants with qualified class\n", "names (``linear.Sersic``, ``operated.Gaussian``); plain names resolve to the standard profiles.\n", diff --git a/notebooks/cluster/lenstool/README.md b/notebooks/cluster/lenstool/README.md index 6d5c5193c..cd1987c82 100644 --- a/notebooks/cluster/lenstool/README.md +++ b/notebooks/cluster/lenstool/README.md @@ -43,7 +43,7 @@ The worked example is **SMACS J0723.3−7327**, the first JWST cluster, using th (median RMS 0.07″, all below 0.29″ — consistent with the published image-plane RMS of 0.32″ through typical magnifications), and the optional Verification II forward-solve checks the image-plane RMS directly. One convention this exercise pinned down the honest way: PyAutoLens's - multi-plane tracer normalizes profile deflections to its **final** plane, so `from_lenstool` + multi-plane tracer normalizes profile deflections to its **final** plane, so `dPIEMass` must be given the highest source redshift in the system as `redshift_source`. ## What's different (and why you might care) diff --git a/notebooks/cluster/lenstool/data.ipynb b/notebooks/cluster/lenstool/data.ipynb index 4edcd1ff7..7ba4322a7 100644 --- a/notebooks/cluster/lenstool/data.ipynb +++ b/notebooks/cluster/lenstool/data.ipynb @@ -26,7 +26,7 @@ " - ``galcat.cat`` \u2192 ``members.csv`` (member catalogue: centres + shape/mag properties,\n", " the ``al.galaxy_table_from_csv`` schema)\n", " - ``best.par`` \u2192 ``mass.csv`` (every optimized ``potential`` section as one\n", - " ``dPIEMassLenstool`` row of the canonical named-galaxy model CSV \u2014 **the .par file as a\n", + " ``dPIEMass`` row of the canonical named-galaxy model CSV \u2014 **the .par file as a\n", " table**, read back with ``al.galaxy_models_from_csv`` like every other cluster dataset)\n", "\n", "__Attribution__\n", @@ -230,7 +230,7 @@ "scaling relation) and 144 cluster members whose parameters Lenstool derived from the scaling\n", "relation. Every section carries the *same five numbers you would type into PyAutoLens*:\n", "\n", - " Lenstool ``potential`` \u2192 ``dPIEMass.from_lenstool`` argument\n", + " Lenstool ``potential`` \u2192 ``dPIEMass`` constructor argument\n", " ---------------------------------------------------------------------\n", " x_centre / y_centre [arcsec] \u2192 centre=(y, x)\n", " ellipticity (a\u00b2-b\u00b2)/(a\u00b2+b\u00b2) \u2192 ellipticity\n", @@ -391,7 +391,7 @@ ")\n", "\n", "# The whole optimized model \u2014 5 named halos + 144 scaling members \u2014 becomes ONE canonical\n", - "# ``mass.csv``: each ``potential`` section is a ``dPIEMassLenstool`` row whose columns are the\n", + "# ``mass.csv``: each ``potential`` section is a ``dPIEMass`` row whose columns are the\n", "# ``.par`` keywords verbatim (sigma, r_core, r_cut, ellipticity, angle_pos) plus the run's\n", "# redshifts and cosmology as flat values. ``modeling.py`` reads it back with the same\n", "# ``al.galaxy_models_from_csv`` call used throughout ``scripts/cluster/``.\n", @@ -404,7 +404,7 @@ "for h in halos:\n", " name = h[\"label\"] if h[\"label\"].startswith(\"O\") else f\"member_{h['label']}\"\n", " profiles_by_galaxy[name] = {\n", - " \"mass\": al.mp.dPIEMassLenstool(\n", + " \"mass\": al.mp.dPIEMass(\n", " centre=(h[\"y\"], h[\"x\"]),\n", " ellipticity=h[\"ellipticity\"],\n", " angle_pos=h[\"angle_pos\"],\n", @@ -425,7 +425,7 @@ " redshifts={name: 0.39 for name in profiles_by_galaxy},\n", ")\n", "\n", - "print(\"Wrote point_datasets.csv, members.csv, mass.csv (149 dPIEMassLenstool rows).\")" + "print(\"Wrote point_datasets.csv, members.csv, mass.csv (149 dPIEMass rows).\")" ], "outputs": [], "execution_count": null diff --git a/notebooks/cluster/lenstool/modeling.ipynb b/notebooks/cluster/lenstool/modeling.ipynb index 6c7993568..dcc119837 100644 --- a/notebooks/cluster/lenstool/modeling.ipynb +++ b/notebooks/cluster/lenstool/modeling.ipynb @@ -44,7 +44,7 @@ "\n", " - **sigma is the fiducial velocity dispersion sigma_LT**, not the physical central velocity\n", " dispersion: sigma_0 = sqrt(3/2) * sigma_LT (Eliasdottir et al. 2007, App. A). PyAutoLens's\n", - " ``from_lenstool`` / ``dPIEMassLenstool`` take sigma_LT \u2014 quote ``v_disp`` values unchanged.\n", + " ``dPIEMass`` takes sigma_LT \u2014 quote ``v_disp`` values unchanged.\n", " Feeding a *measured* stellar velocity dispersion here overestimates the mass by 50%.\n", " - **The x axis points West** in Lenstool's relative frame (data.py verifies this against the\n", " data). All coordinates in this script live in that frame, so numbers compare directly to the\n", @@ -62,7 +62,7 @@ "__Contents__\n", "\n", "- **Load Data:** the CSVs written by ``data.py``.\n", - "- **The Published Model, Reconstructed:** 149 ``from_lenstool`` profiles + 21 point sources.\n", + "- **The Published Model, Reconstructed:** 149 Lenstool-native ``dPIEMass`` profiles + 21 point sources.\n", "- **Verification I \u2014 source-plane compactness:** observed images trace to tight source groups.\n", "- **Verification II \u2014 image-plane RMS (optional):** forward-solve vs the published 0.32\".\n", "- **Critical Curves (optional):** the per-source-plane critical curves over the HST image.\n", @@ -169,7 +169,7 @@ " exist, the model-optimized values of ``best.par`` otherwise).\n", " - ``mass.csv`` \u2014 the complete optimized mass model in the **canonical named-galaxy CSV**\n", " (the same ``al.galaxy_models_from_csv`` format every cluster script uses): 149 rows of\n", - " ``profile_class = dPIEMassLenstool``, one per ``potential`` section of ``best.par``, whose\n", + " ``profile_class = dPIEMass``, one per ``potential`` section of ``best.par``, whose\n", " columns are the ``.par`` keywords verbatim. The five individually-optimized halos are named\n", " O1 (cluster-scale), O2 (BCG), O3 (\"dNW\"), O4 (\"ICL\"), O5 (\"eCM\"); the 144 scaling members\n", " are ``member_``.\n", @@ -193,7 +193,7 @@ "members_table = al.galaxy_table_from_csv(file_path=dataset_path / \"members.csv\")\n", "\n", "print(\n", - " f\"mass.csv: {len(mass_table.rows)} dPIEMassLenstool rows | \"\n", + " f\"mass.csv: {len(mass_table.rows)} dPIEMass rows | \"\n", " f\"members.csv: {len(members_table.luminosities)} catalogue members\"\n", ")" ], @@ -206,7 +206,8 @@ "source": [ "__The Published Model, Reconstructed__\n", "\n", - "Every ``potential`` section of ``best.par`` becomes one ``al.mp.dPIEMass`` via ``from_lenstool``\n", + "Every ``potential`` section of ``best.par`` becomes one ``al.mp.dPIEMass`` \u2014 the default dPIE is\n", + "Lenstool's native parameterization, so the ``.par`` keywords are the constructor arguments\n", "\u2014 the arguments are the ``.par`` keywords, verbatim. This is the whole point of the Lenstool-native\n", "API: nothing is transcribed by hand, and the sqrt(3/2) sigma convention, the ellipticity\n", "conversion and the D_LS/D_S normalization are handled (and unit-tested) inside PyAutoLens.\n", @@ -227,7 +228,7 @@ "source": [ "Z_REF_SOURCE = max(float(dataset.redshift) for dataset in dataset_list)\n", "\n", - "# One call: every mass.csv row instantiates its dPIEMassLenstool with the .par values \u2014\n", + "# One call: every mass.csv row instantiates its dPIEMass with the .par values \u2014\n", "# the redshift_source (final-plane) normalization and the run's H0/Om0 travel inside the\n", "# CSV columns, so nothing here needs to remember them.\n", "lens_galaxies = list(al.galaxies_from_csv_tables(mass_table).values())\n", @@ -402,7 +403,7 @@ "\n", "Everything above used the published answer. A real analysis *fits*: the model below reproduces the\n", "composition Lenstool optimized (``input.par``), using the Lenstool-parameterized profile\n", - "``al.mp.dPIEMassLenstool`` so every free parameter, prior bound and posterior number is in\n", + "``al.mp.dPIEMass`` so every free parameter, prior bound and posterior number is in\n", "Lenstool units:\n", "\n", " - **O1, cluster halo**: centre U(-5,5)\" (both axes), ellipticity U(0,0.8), angle U(-90,90),\n", @@ -415,7 +416,7 @@ " r_cut U(50,1000) kpc, sigma U(0,700). [7 free]\n", " - **O5 \"eCM\"**: centre fixed; ellipticity U(0,0.6), angle U(-90,90), r_core U(0,10) kpc,\n", " r_cut U(10,200) kpc, sigma U(0.1,300). [5 free]\n", - " - **potfile members**: every catalogue member gets a ``dPIEMassLenstool`` with centre, shape and\n", + " - **potfile members**: every catalogue member gets a ``dPIEMass`` with centre, shape and\n", " angle *fixed to the light* (``galcat.cat``) and its sigma / r_cut derived from two shared free\n", " parameters exactly as ``potfile`` defines \u2014\n", "\n", @@ -453,7 +454,7 @@ " members_table.properties[\"ellipticity\"],\n", " members_table.properties[\"angle_pos\"],\n", "):\n", - " mass = af.Model(al.mp.dPIEMassLenstool)\n", + " mass = af.Model(al.mp.dPIEMass)\n", " mass.centre = tuple(centre)\n", " mass.ellipticity = ellipticity\n", " mass.angle_pos = angle_pos\n", diff --git a/notebooks/cluster/likelihood_function.ipynb b/notebooks/cluster/likelihood_function.ipynb index f2cdc8322..2e9eb6d3a 100644 --- a/notebooks/cluster/likelihood_function.ipynb +++ b/notebooks/cluster/likelihood_function.ipynb @@ -225,22 +225,28 @@ "# REFERENCE_LUMINOSITY is an explicit fixed constant (Lenstool's \"mag0\"), not the\n", "# sample max, and matches the simulator truth so members are reproduced exactly.\n", "scaling_galaxies = []\n", - "SCALING_B0_REF_TRUTH = 0.190\n", - "SCALING_EXPONENT = 0.5\n", - "SCALING_RA_REF = 0.158\n", - "SCALING_RS_REF = 15.8\n", + "SCALING_SIGMA_REF_TRUTH = 85.0\n", + "SCALING_SIGMA_EXPONENT = 0.25\n", + "SCALING_RADIUS_EXPONENT = 0.5\n", + "SCALING_R_CORE_REF = 0.158\n", + "SCALING_R_CUT_REF = 15.8\n", "REFERENCE_LUMINOSITY = 1.0\n", "for centre, luminosity in zip(\n", " scaling_table.centres.in_list, scaling_table.luminosities\n", "):\n", " luminosity_ratio = luminosity / REFERENCE_LUMINOSITY\n", - " ra = SCALING_RA_REF * luminosity_ratio**SCALING_EXPONENT\n", - " b0 = SCALING_B0_REF_TRUTH * luminosity_ratio**SCALING_EXPONENT\n", - " rs = SCALING_RS_REF * luminosity_ratio**SCALING_EXPONENT\n", " scaling_galaxies.append(\n", " al.Galaxy(\n", " redshift=redshift_lens,\n", - " mass=al.mp.dPIEMassSph(centre=tuple(centre), ra=ra, rs=rs, b0=b0),\n", + " mass=al.mp.dPIEMassSph(\n", + " centre=tuple(centre),\n", + " sigma=SCALING_SIGMA_REF_TRUTH\n", + " * luminosity_ratio**SCALING_SIGMA_EXPONENT,\n", + " r_core=SCALING_R_CORE_REF * luminosity_ratio**SCALING_RADIUS_EXPONENT,\n", + " r_cut=SCALING_R_CUT_REF * luminosity_ratio**SCALING_RADIUS_EXPONENT,\n", + " redshift_object=redshift_lens,\n", + " redshift_source=max(source_redshifts),\n", + " ),\n", " )\n", " )" ], @@ -256,8 +262,9 @@ "The tracer carries:\n", "\n", " - 2 main lens galaxies (BCG + satellite) \u2014 individually-modelled dPIE mass profiles.\n", - " - 10 scaling-tier member galaxies \u2014 dPIE mass profiles whose ``ra``, ``rs`` and ``b0`` derive from the\n", - " reference-anchored scaling relation ``b0 = b0_ref \u00d7 (L/L_ref)^0.5`` (Lenstool convention).\n", + " - 10 scaling-tier member galaxies \u2014 dPIE mass profiles whose ``sigma``, ``r_core`` and ``r_cut`` derive\n", + " from the reference-anchored scaling relation ``sigma = sigma_ref \u00d7 (L/L_ref)^0.25``,\n", + " ``radii \u221d (L/L_ref)^0.5`` (Lenstool convention).\n", " - 1 host dark matter halo \u2014 ``NFWMCRLudlowSph`` at the cluster centre.\n", " - 2 source galaxies \u2014 ``Point`` profiles at distinct redshifts (multi-plane).\n", "\n", diff --git a/notebooks/cluster/modeling.ipynb b/notebooks/cluster/modeling.ipynb index d0f548eed..3edba473b 100644 --- a/notebooks/cluster/modeling.ipynb +++ b/notebooks/cluster/modeling.ipynb @@ -48,8 +48,8 @@ "\n", " - There are 2 main lens galaxies with ``dPIEMassSph`` total mass distributions, each with their centre\n", " fixed to the values written out by the simulator [6 parameters].\n", - " - There are 10 scaling-tier member galaxies. Each carries a ``dPIEMassSph`` mass with centre fixed,\n", - " ``ra`` fixed, and ``b0`` / ``rs`` derived from the shared reference-anchored scaling relation and the\n", + " - There are 10 scaling-tier member galaxies. Each carries a ``dPIEMassSph`` mass with centre fixed and\n", + " ``sigma`` / ``r_core`` / ``r_cut`` derived from the shared reference-anchored scaling relation and the\n", " per-member luminosity [1 parameter total for the entire tier].\n", " - There is 1 standalone ``NFWMCRLudlowSph`` host dark matter halo with its centre fixed and a free\n", " ``mass_at_200`` [1 parameter].\n", @@ -58,11 +58,16 @@ "\n", "The number of free parameters and therefore the dimensionality of non-linear parameter space is N=12.\n", "\n", + "The ``dPIEMassSph`` profile is parameterized in Lenstool's native convention \u2014 free parameters\n", + "``sigma`` (fiducial velocity dispersion ``v_disp`` in km/s), ``r_core`` and ``r_cut`` (arcsec) \u2014 so the\n", + "fitted posterior reads like a Lenstool results table (see ``cluster/simulator.py`` __dPIE Mass Profile__\n", + "for the conventions, and ``cluster/lenstool/`` for a published-model worked example).\n", + "\n", "The defining feature of cluster modeling is the scaling tier: 10 lower-mass members are fit jointly\n", - "with a *single free parameter* (``b0_ref``, the lens strength of a galaxy at the reference magnitude;\n", - "the relation's exponent is fixed at the Faber-Jackson value of 0.5, following the Lenstool convention).\n", - "Adding more members to ``scaling_galaxies.csv`` in the future does not grow the dimensionality of\n", - "parameter space.\n", + "with a *single free parameter* (``sigma_ref``, the fiducial velocity dispersion of a galaxy at the\n", + "reference magnitude; the relation's exponents are fixed at the Faber-Jackson values, following the\n", + "Lenstool convention). Adding more members to ``scaling_galaxies.csv`` in the future does not grow the\n", + "dimensionality of parameter space.\n", "\n", "__Simulation__\n", "\n", @@ -386,10 +391,10 @@ " with a ``dPIEMassSph`` total mass profile whose centre is fixed to ``main_lens_centres[i]``.\n", "\n", " - ``scaling_galaxies``: The 10 scaling-tier cluster members. Each carries a ``dPIEMassSph`` mass with\n", - " centre fixed (from the CSV), ``ra`` fixed, and ``b0`` / ``rs`` derived from the shared\n", - " reference-anchored scaling relation (single free normalization ``b0_ref``, exponents fixed at 0.5)\n", - " plus the per-member luminosity. The whole tier contributes 1 free parameter to the model regardless\n", - " of how many members are in the CSV.\n", + " centre fixed (from the CSV) and ``sigma`` / ``r_core`` / ``r_cut`` derived from the shared\n", + " reference-anchored scaling relation (single free normalization ``sigma_ref``; exponents fixed at\n", + " 0.25 for sigma and 0.5 for the radii) plus the per-member luminosity. The whole tier contributes\n", + " 1 free parameter to the model regardless of how many members are in the CSV.\n", "\n", " - ``host_halo``: A single standalone ``Galaxy`` carrying the cluster's ``NFWMCRLudlowSph`` dark matter\n", " halo. The halo is *not* tied to any individual member \u2014 it sits \"on top of\" the members and\n", @@ -423,12 +428,12 @@ "We compose a lens model where:\n", "\n", " - The 2 main lens galaxies each have a ``dPIEMassSph`` mass profile with centre fixed and free\n", - " ``ra``, ``rs``, ``b0`` \u2014 3 free parameters per galaxy [6 parameters].\n", - " - The 10 scaling-tier members share a single free parameter: ``b0_ref``, the lens strength of a\n", - " galaxy at the reference magnitude. Each member's ``ra``, ``rs`` and ``b0`` are computed as\n", - " ``ra_ref * (L / L_ref) ** 0.5``, ``rs_ref * (L / L_ref) ** 0.5`` and ``b0_ref * (L / L_ref) ** 0.5``\n", - " with the exponent fixed; ``ra_ref`` (0.158\") and ``rs_ref`` (15.8\") are held fixed at the simulator\n", - " truth values [1 parameter].\n", + " ``sigma``, ``r_core``, ``r_cut`` \u2014 3 free parameters per galaxy [6 parameters].\n", + " - The 10 scaling-tier members share a single free parameter: ``sigma_ref``, the fiducial velocity\n", + " dispersion of a galaxy at the reference magnitude. Each member's parameters are computed as\n", + " ``sigma_ref * (L / L_ref) ** 0.25``, ``r_core_ref * (L / L_ref) ** 0.5`` and\n", + " ``r_cut_ref * (L / L_ref) ** 0.5`` with the exponents fixed; ``r_core_ref`` (0.158\") and\n", + " ``r_cut_ref`` (15.8\") are held fixed at the simulator truth values [1 parameter].\n", " - The host halo has an ``NFWMCRLudlowSph`` mass profile with centre fixed and a free ``mass_at_200``\n", " [1 parameter].\n", " - Each source has a ``Point`` model with free ``centre_0`` / ``centre_1`` priors initialised from the\n", @@ -442,30 +447,31 @@ "published cluster strong-lensing analysis (Limousin et al. 2005; Eliasdottir et al. 2007; Bergamini et\n", "al. 2019):\n", "\n", - " ra_i = ra_ref * (L_i / L_ref) ** 0.5\n", - " rs_i = rs_ref * (L_i / L_ref) ** 0.5\n", - " b0_i = b0_ref * (L_i / L_ref) ** 0.5\n", + " sigma_i = sigma_ref * (L_i / L_ref) ** 0.25\n", + " r_core_i = r_core_ref * (L_i / L_ref) ** 0.5\n", + " r_cut_i = r_cut_ref * (L_i / L_ref) ** 0.5\n", "\n", "The reference luminosity ``L_ref`` is an **explicit fixed constant** (Lenstool's reference magnitude\n", "``mag0``), *not* the maximum luminosity of the current sample. Anchoring to a fixed reference makes the\n", "normalization invariant to which galaxies are placed in the scaling tier, and gives the single free\n", - "parameter ``b0_ref`` a stable, interpretable meaning: the lens strength (roughly the Einstein radius) of\n", + "parameter ``sigma_ref`` a stable, interpretable meaning: the fiducial velocity dispersion (in km/s) of\n", "a galaxy *at the reference magnitude*, for which a prior range is easy to motivate \u2014 unlike an abstract\n", "multiplicative factor whose units depend on the (arbitrary) luminosity normalization. In a real analysis\n", "set ``L_ref`` to the BCG magnitude (or a characteristic L*); here we use a fiducial ``L_ref = 1.0``. Only\n", "luminosity *ratios* enter, so the CSV's luminosity units are irrelevant; magnitude catalogues convert via\n", "``L_i / L_ref = 10 ** (0.4 * (m_ref - m_i))``.\n", "\n", - "The exponent is *fixed* at 0.5 rather than fitted: b0 \u221d sigma\u00b2 for the dPIE, and Faber-Jackson\n", - "(sigma \u221d L^(1/4)) gives b0 \u221d L^(1/2). The core radius ``ra`` and truncation radius ``rs`` scale with the\n", - "same fixed exponent, mirroring Lenstool's r_core \u221d L^(1/2) and r_cut \u221d L^(1/2). Freeing the exponent (or\n", - "``rs_ref`` / ``ra_ref``) is a one-line change shown in the code comment below \u2014 useful as a systematics\n", - "test, at the cost of the degeneracy between normalization and slope that the fixed-exponent convention\n", - "exists to avoid. When member velocity dispersions are available, the standard refinement is to calibrate\n", - "the exponents kinematically (Bergamini et al. 2019: sigma \u221d L^0.27-0.28 from MUSE member kinematics, i.e.\n", - "a b0 exponent \u2248 0.55, with the r_cut exponent from the fundamental plane).\n", - "\n", - "The simulator's truth value is ``b0_ref = 0.190`` arcsec (at ``L_ref = 1.0``). The prior below is much\n", + "The exponents are *fixed* rather than fitted: Faber-Jackson (L \u221d sigma^4) gives sigma \u221d L^(1/4), and\n", + "constant mass-to-light (M \u221d sigma^2 * r_cut \u221d L) gives r_cut \u221d L^(1/2) \u2014 exactly Lenstool's ``potfile``\n", + "scalings (``vdslope 4``, ``slope 4``), applied to the core radius too. Since the dPIE lens strength obeys\n", + "b0 \u221d sigma^2, this is equivalent to the b0 \u221d L^(1/2) scaling of the internal parameterization. Freeing an\n", + "exponent (or ``r_cut_ref`` / ``r_core_ref``) is a one-line change shown in the code comment below \u2014\n", + "useful as a systematics test, at the cost of the degeneracy between normalization and slope that the\n", + "fixed-exponent convention exists to avoid. When member velocity dispersions are available, the standard\n", + "refinement is to calibrate the exponents kinematically (Bergamini et al. 2019: sigma \u221d L^0.27-0.28 from\n", + "MUSE member kinematics, with the r_cut exponent from the fundamental plane).\n", + "\n", + "The simulator's truth value is ``sigma_ref = 85.0`` km/s (at ``L_ref = 1.0``). The prior below is much\n", "wider than the truth to give the search room." ] }, @@ -483,11 +489,16 @@ "\n", "galaxy_models = al.galaxy_af_models_from_csv_tables(mass_table, point_table)\n", "\n", - "# Main Lens Galaxies: free dPIE ra / rs / b0 on each; centre stays fixed at the CSV value.\n", + "# Main Lens Galaxies: free dPIE sigma / r_core / r_cut on each; centre and redshifts\n", + "# stay fixed at the CSV values, and the cosmology constants H0 / Om0 are pinned (they\n", + "# are model *constants*, not parameters to sample \u2014 if left unset they would inherit\n", + "# the config's default priors and float).\n", "for name in (\"lens_0\", \"lens_1\"):\n", - " galaxy_models[name].mass.ra = af.UniformPrior(lower_limit=1.0, upper_limit=15.0)\n", - " galaxy_models[name].mass.rs = af.UniformPrior(lower_limit=5.0, upper_limit=40.0)\n", - " galaxy_models[name].mass.b0 = af.UniformPrior(lower_limit=0.1, upper_limit=10.0)\n", + " galaxy_models[name].mass.sigma = af.UniformPrior(lower_limit=50.0, upper_limit=600.0)\n", + " galaxy_models[name].mass.r_core = af.UniformPrior(lower_limit=1.0, upper_limit=15.0)\n", + " galaxy_models[name].mass.r_cut = af.UniformPrior(lower_limit=5.0, upper_limit=40.0)\n", + " galaxy_models[name].mass.H0 = 67.66\n", + " galaxy_models[name].mass.Om0 = 0.30966\n", "\n", "# Host Halo: free mass_at_200; centre + redshift_object + redshift_source stay fixed.\n", "galaxy_models[\"host_halo\"].dark.mass_at_200 = af.LogUniformPrior(\n", @@ -509,7 +520,7 @@ " mean=float(np.mean(positions[:, 1])), sigma=3.0\n", " )\n", "\n", - "# Scaling Tier Members (dPIEMassSph; ra, rs and b0 all derived from the\n", + "# Scaling Tier Members (dPIEMassSph; sigma, r_core and r_cut all derived from the\n", "# reference-anchored scaling relation).\n", "#\n", "# The reference luminosity is an EXPLICIT FIXED constant (Lenstool's reference\n", @@ -518,23 +529,25 @@ "# magnitude in a real analysis; here it is a fiducial L* = 1.0.\n", "reference_luminosity = 1.0\n", "\n", - "# b0_ref is defined ONCE outside the loop \u2014 the tier's only free parameter, the\n", - "# lens strength of a reference-magnitude galaxy. Every member's ra, rs and b0 are\n", - "# derived by scaling the fixed reference values (ra_ref, rs_ref) or the free b0_ref\n", - "# by its luminosity ratio, with the exponent fixed at the Faber-Jackson value 0.5.\n", - "# The entire tier therefore contributes 1 free parameter regardless of how many\n", + "# sigma_ref is defined ONCE outside the loop \u2014 the tier's only free parameter, the\n", + "# fiducial velocity dispersion of a reference-magnitude galaxy (km/s, so the prior\n", + "# range below is physically interpretable). Every member's sigma, r_core and r_cut\n", + "# are derived by scaling the free sigma_ref or the fixed reference radii by its\n", + "# luminosity ratio, with the exponents fixed at the Faber-Jackson values\n", + "# (sigma \u221d L^0.25, radii \u221d L^0.5 \u2014 Lenstool's potfile vdslope/slope = 4). The\n", + "# entire tier therefore contributes 1 free parameter regardless of how many\n", "# members are in scaling_galaxies.csv.\n", "#\n", - "# To free the exponent as a systematics test, replace the fixed value with e.g.\n", - "# `scaling_exponent = af.UniformPrior(lower_limit=0.0, upper_limit=1.0)` \u2014 every\n", - "# member's b0 then derives from two shared parameters, as in older versions of\n", - "# this example.\n", + "# To free an exponent as a systematics test, replace the fixed value with e.g.\n", + "# `scaling_sigma_exponent = af.UniformPrior(lower_limit=0.0, upper_limit=0.5)` \u2014\n", + "# every member's sigma then derives from two shared parameters.\n", "\n", - "scaling_b0_ref = af.UniformPrior(lower_limit=0.0, upper_limit=1.0)\n", - "scaling_exponent = 0.5\n", + "scaling_sigma_ref = af.UniformPrior(lower_limit=0.0, upper_limit=200.0)\n", + "scaling_sigma_exponent = 0.25\n", + "scaling_radius_exponent = 0.5\n", "\n", - "scaling_ra_ref_fixed = 0.158\n", - "scaling_rs_ref_fixed = 15.8\n", + "scaling_r_core_ref_fixed = 0.158\n", + "scaling_r_cut_ref_fixed = 15.8\n", "\n", "scaling_galaxies_list = []\n", "for centre, luminosity in zip(\n", @@ -544,9 +557,13 @@ "\n", " mass = af.Model(al.mp.dPIEMassSph)\n", " mass.centre = tuple(centre)\n", - " mass.ra = scaling_ra_ref_fixed * luminosity_ratio**scaling_exponent\n", - " mass.rs = scaling_rs_ref_fixed * luminosity_ratio**scaling_exponent\n", - " mass.b0 = scaling_b0_ref * luminosity_ratio**scaling_exponent\n", + " mass.sigma = scaling_sigma_ref * luminosity_ratio**scaling_sigma_exponent\n", + " mass.r_core = scaling_r_core_ref_fixed * luminosity_ratio**scaling_radius_exponent\n", + " mass.r_cut = scaling_r_cut_ref_fixed * luminosity_ratio**scaling_radius_exponent\n", + " mass.redshift_object = redshift_lens\n", + " mass.redshift_source = max(source_redshifts)\n", + " mass.H0 = 67.66\n", + " mass.Om0 = 0.30966\n", "\n", " scaling_galaxies_list.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass))\n", "\n", diff --git a/notebooks/cluster/simulator.ipynb b/notebooks/cluster/simulator.ipynb index 26f91cf16..24e7df1cc 100644 --- a/notebooks/cluster/simulator.ipynb +++ b/notebooks/cluster/simulator.ipynb @@ -72,9 +72,9 @@ " relation. Each member is individually weak compared to the main galaxies or the host halo, but the\n", " population together perturbs the deflection field non-trivially \u2014 exactly the regime in which the\n", " scaling-relation tier of the modeling API earns its keep. The number of free parameters does not grow\n", - " with the number of scaling members; a single shared normalization `b0_ref` (the lens strength of a\n", - " galaxy at the reference magnitude, with the relation's exponent fixed at the Faber-Jackson value)\n", - " determines every member's mass from its luminosity.\n", + " with the number of scaling members; a single shared normalization `sigma_ref` (the fiducial velocity\n", + " dispersion of a galaxy at the reference magnitude, with the relation's exponents fixed at the\n", + " Faber-Jackson values) determines every member's mass from its luminosity.\n", "\n", "- `host_halo_galaxy`: A standalone `Galaxy` holding the cluster's `NFWMCRLudlowSph` dark matter halo. It\n", " is not tied to any individual member galaxy \u2014 the halo is a separate mass component sitting \"on top of\"\n", @@ -92,49 +92,62 @@ "\n", "The cluster member galaxies use the dual Pseudo-Isothermal Elliptical (dPIE) mass profile introduced in\n", "Eliasdottir 2007 (https://arxiv.org/abs/0710.5636), the de facto standard for cluster strong lens modeling.\n", - "In spherical form (`dPIEMassSph`), its parameters are:\n", - "\n", - " - `ra` (arcsec): the core radius, below which the density profile flattens (kept small, ~0.05\u20130.1\" at z=0.5).\n", - " - `rs` (arcsec): the truncation radius, above which the density falls as R^-4 (kept ~10\u201330\" for cluster members).\n", - " - `b0` (arcsec): the mass normalization, roughly setting the galaxy-scale Einstein radius.\n", + "In spherical form (`dPIEMassSph`), its parameters are **Lenstool's native ones** \u2014 the same numbers that\n", + "appear in a Lenstool ``.par`` file or a published cluster paper's results table:\n", + "\n", + " - `sigma` (km/s): the fiducial velocity dispersion (Lenstool's ``v_disp``, sigma_LT). This is *not* the\n", + " physical central dispersion \u2014 sigma_0 = sqrt(3/2) * sigma_LT (Eliasdottir et al. 2007, App. A) \u2014 so\n", + " quote ``.par``/table values unchanged and never feed a measured stellar dispersion here.\n", + " - `r_core` (arcsec): the core radius, below which the density profile flattens.\n", + " - `r_cut` (arcsec): the truncation radius, above which the density falls as R^-4 (~10\u201330\" for cluster members).\n", + " - `redshift_object` / `redshift_source`: the redshifts entering the D_LS / D_S normalization of the\n", + " internal lens strength ``b0 = 6 * 648000 * (sigma / c)^2 * (D_LS / D_S)`` arcsec. In this multi-plane\n", + " cluster, ``redshift_source`` is anchored to the *furthest* source (``z = 2.0``) \u2014 the same convention\n", + " as the ``NFWMCRLudlowSph`` host halo below, because the tracer normalizes profile deflections to the\n", + " final source plane.\n", + "\n", + "The internal ``(ra, rs, b0)`` parameterization remains available as ``dPIEMassB0Sph`` (non-standard; b0\n", + "is the lens strength in arcsec, with b0 proportional to sigma^2).\n", "\n", "Per-galaxy values for the 2 main-tier galaxies are hand-tuned below; for the 10 scaling-tier members they\n", "are derived from each member's luminosity via the relation described next.\n", "\n", "__Luminosity-Mass Scaling Relation__\n", "\n", - "The 10 scaling members share a reference-anchored relation for the dPIE mass normalization \u2014 the\n", - "convention used by Lenstool and essentially every published cluster strong-lensing analysis\n", - "(Limousin et al. 2005; Eliasdottir et al. 2007; Bergamini et al. 2019):\n", + "The 10 scaling members share a reference-anchored relation \u2014 the convention used by Lenstool and\n", + "essentially every published cluster strong-lensing analysis (Limousin et al. 2005; Eliasdottir et al.\n", + "2007; Bergamini et al. 2019):\n", "\n", - " ra_i = ra_ref * (L_i / L_ref) ** 0.5\n", - " rs_i = rs_ref * (L_i / L_ref) ** 0.5\n", - " b0_i = b0_ref * (L_i / L_ref) ** 0.5\n", + " sigma_i = sigma_ref * (L_i / L_ref) ** 0.25\n", + " r_core_i = r_core_ref * (L_i / L_ref) ** 0.5\n", + " r_cut_i = r_cut_ref * (L_i / L_ref) ** 0.5\n", "\n", "where ``L_ref`` is an **explicit fixed reference luminosity** (Lenstool's reference magnitude ``mag0``),\n", - "*not* the maximum luminosity of the sample, and ``b0_ref`` is the lens strength of a galaxy *at that\n", - "reference magnitude*. Anchoring to a fixed reference makes the normalization physically interpretable \u2014\n", - "the Einstein-radius-like strength of a reference-magnitude galaxy \u2014 and invariant to which galaxies are\n", - "placed in the tier, which is what makes a sensible prior range easy to define. The exponent is **fixed\n", - "at 0.5** rather than fitted: for the dPIE, ``b0`` is proportional to the velocity dispersion squared, and\n", - "the Faber-Jackson relation (L \u221d sigma^4, i.e. sigma \u221d L^(1/4)) then gives b0 \u221d L^(1/2). Lenstool applies\n", - "the same fixed-exponent scaling to the core and truncation radii (r_core \u221d L^(1/2), r_cut \u221d L^(1/2)),\n", - "which is why ``ra`` and ``rs`` scale here too.\n", - "\n", - "Truth values used in this simulator are ``b0_ref = 0.190`` arcsec, ``rs_ref = 15.8`` arcsec and\n", - "``ra_ref = 0.158`` arcsec at a fiducial reference luminosity ``L_ref = 1.0`` (the model in\n", - "``modeling.py`` / ``start_here.py`` uses the SAME reference, so members stay consistent by construction). Luminosities are log-spaced across roughly 0.05\u20130.40, so\n", - "per-member ``b0`` values run from ~0.042 to 0.12 arcsec \u2014 each member is individually well below the\n", - "BCG (``b0 = 3.0``) but the 10 of them together perturb the deflection field by ~10\u201315%.\n", - "\n", - "The modeling script promotes ``b0_ref`` to the tier's single free parameter and recovers the truth value\n", - "when fit to the simulated point datasets. Adding more scaling members amounts to adding rows to\n", + "*not* the maximum luminosity of the sample, and ``sigma_ref`` is the fiducial velocity dispersion of a\n", + "galaxy *at that reference magnitude*. Anchoring to a fixed reference makes the normalization physically\n", + "interpretable \u2014 the velocity dispersion of a reference-magnitude galaxy, in km/s \u2014 and invariant to which\n", + "galaxies are placed in the tier, which is what makes a sensible prior range easy to define. The exponents\n", + "are **fixed** rather than fitted: the Faber-Jackson relation (L \u221d sigma^4) gives sigma \u221d L^(1/4), and\n", + "constant mass-to-light (M \u221d sigma^2 * r_cut \u221d L) then gives r_cut \u221d L^(1/2) \u2014 Lenstool's ``potfile``\n", + "applies exactly these scalings (``vdslope 4``, ``slope 4``), including to the core radius. Because the\n", + "lens strength b0 is proportional to sigma^2, this is equivalent to the b0 \u221d L^(1/2) scaling of the\n", + "internal parameterization.\n", + "\n", + "Truth values used in this simulator are ``sigma_ref = 85.0`` km/s, ``r_cut_ref = 15.8`` arcsec and\n", + "``r_core_ref = 0.158`` arcsec at a fiducial reference luminosity ``L_ref = 1.0`` (the model in\n", + "``modeling.py`` / ``start_here.py`` uses the SAME reference, so members stay consistent by construction).\n", + "Luminosities are log-spaced across roughly 0.05\u20130.40, so per-member dispersions run from ~40 to ~68 km/s\n", + "\u2014 each member is individually well below the BCG (``sigma = 330`` km/s) but the 10 of them together\n", + "perturb the deflection field by ~10\u201315%.\n", + "\n", + "The modeling script promotes ``sigma_ref`` to the tier's single free parameter and recovers the truth\n", + "value when fit to the simulated point datasets. Adding more scaling members amounts to adding rows to\n", "``scaling_galaxies.csv`` \u2014 the tier's free-parameter count stays at 1. Note that only the luminosity\n", "*ratios* ``L_i / L_ref`` enter the relation, so the units of the luminosity column are irrelevant;\n", "observational catalogues quoting magnitudes convert via ``L_i / L_ref = 10 ** (0.4 * (m_ref - m_i))``.\n", "Kinematic calibrations of the exponent exist for when higher fidelity is needed \u2014 Bergamini et al. 2019\n", - "measure sigma \u221d L^0.27-0.28 from MUSE member kinematics (b0 exponent \u2248 0.55) and derive the r_cut\n", - "exponent from the fundamental plane \u2014 but 0.5 is the standard default.\n", + "measure sigma \u221d L^0.27-0.28 from MUSE member kinematics and derive the r_cut exponent from the\n", + "fundamental plane (beta_cut = gamma - 2*alpha + 1) \u2014 but 0.25 / 0.5 are the standard defaults.\n", "\n", "__NFWMCRLudlow Host Halo__\n", "\n", @@ -409,9 +422,9 @@ "metadata": {}, "source": [ "main_lens_dpie_params = [\n", - " # (ra, rs, b0) per galaxy \u2014 arcsec\n", - " (8.0, 20.0, 3.0), # BCG \u2014 strongest\n", - " (5.0, 12.0, 1.2), # satellite\n", + " # (r_core, r_cut, sigma) per galaxy \u2014 arcsec, arcsec, km/s\n", + " (8.0, 20.0, 330.0), # BCG \u2014 strongest\n", + " (5.0, 12.0, 210.0), # satellite\n", "]\n", "\n", "main_lens_sersic_params = [\n", @@ -421,7 +434,7 @@ "]\n", "\n", "main_lens_galaxies = []\n", - "for centre, (ra, rs, b0), (intensity, effective_radius, sersic_index) in zip(\n", + "for centre, (r_core, r_cut, sigma), (intensity, effective_radius, sersic_index) in zip(\n", " main_lens_centres, main_lens_dpie_params, main_lens_sersic_params\n", "):\n", " bulge = al.lp.SersicSph(\n", @@ -430,7 +443,14 @@ " effective_radius=effective_radius,\n", " sersic_index=sersic_index,\n", " )\n", - " mass = al.mp.dPIEMassSph(centre=centre, ra=ra, rs=rs, b0=b0)\n", + " mass = al.mp.dPIEMassSph(\n", + " centre=centre,\n", + " sigma=sigma,\n", + " r_core=r_core,\n", + " r_cut=r_cut,\n", + " redshift_object=redshift_lens,\n", + " redshift_source=max(source_redshifts),\n", + " )\n", " main_lens_galaxies.append(al.Galaxy(redshift=redshift_lens, bulge=bulge, mass=mass))" ], "outputs": [], @@ -444,24 +464,25 @@ "\n", "The 10 cluster members modelled collectively via the luminosity-mass scaling relation (see the\n", "``__Luminosity-Mass Scaling Relation__`` section of the module docstring). The simulator hardcodes the\n", - "truth value of ``b0_ref`` (the lens strength of a galaxy at the reference magnitude) and derives each\n", - "member's ``ra``, ``rs`` and ``b0`` from its luminosity ratio to the reference, with the exponent fixed\n", - "at the Faber-Jackson value of 0.5. The reference luminosity is an EXPLICIT FIXED constant (Lenstool's\n", - "``mag0``), NOT the sample max \u2014 here a fiducial ``L_ref = 1.0``; the model in ``modeling.py`` /\n", - "``start_here.py`` uses the SAME reference so members stay consistent by construction. Light profiles use\n", - "the per-member luminosity as the central intensity so the rendered image visibly traces the scaling-tier\n", - "population." + "truth value of ``sigma_ref`` (the fiducial velocity dispersion of a galaxy at the reference magnitude)\n", + "and derives each member's ``sigma``, ``r_core`` and ``r_cut`` from its luminosity ratio to the reference,\n", + "with the exponents fixed at the Faber-Jackson values (sigma \u221d L^0.25; radii \u221d L^0.5). The reference\n", + "luminosity is an EXPLICIT FIXED constant (Lenstool's ``mag0``), NOT the sample max \u2014 here a fiducial\n", + "``L_ref = 1.0``; the model in ``modeling.py`` / ``start_here.py`` uses the SAME reference so members stay\n", + "consistent by construction. Light profiles use the per-member luminosity as the central intensity so the\n", + "rendered image visibly traces the scaling-tier population." ] }, { "cell_type": "code", "metadata": {}, "source": [ - "scaling_b0_ref_truth = 0.190\n", - "scaling_exponent = 0.5\n", + "scaling_sigma_ref_truth = 85.0\n", + "scaling_sigma_exponent = 0.25\n", + "scaling_radius_exponent = 0.5\n", "reference_luminosity = 1.0\n", - "scaling_ra_ref = 0.158\n", - "scaling_rs_ref = 15.8\n", + "scaling_r_core_ref = 0.158\n", + "scaling_r_cut_ref = 15.8\n", "\n", "scaling_galaxies = []\n", "for centre, luminosity in zip(scaling_galaxies_centres, scaling_galaxies_luminosities):\n", @@ -472,10 +493,14 @@ " sersic_index=3.0,\n", " )\n", " luminosity_ratio = luminosity / reference_luminosity\n", - " ra = scaling_ra_ref * luminosity_ratio**scaling_exponent\n", - " b0 = scaling_b0_ref_truth * luminosity_ratio**scaling_exponent\n", - " rs = scaling_rs_ref * luminosity_ratio**scaling_exponent\n", - " mass = al.mp.dPIEMassSph(centre=centre, ra=ra, rs=rs, b0=b0)\n", + " mass = al.mp.dPIEMassSph(\n", + " centre=centre,\n", + " sigma=scaling_sigma_ref_truth * luminosity_ratio**scaling_sigma_exponent,\n", + " r_core=scaling_r_core_ref * luminosity_ratio**scaling_radius_exponent,\n", + " r_cut=scaling_r_cut_ref * luminosity_ratio**scaling_radius_exponent,\n", + " redshift_object=redshift_lens,\n", + " redshift_source=max(source_redshifts),\n", + " )\n", " scaling_galaxies.append(al.Galaxy(redshift=redshift_lens, bulge=bulge, mass=mass))" ], "outputs": [], @@ -751,8 +776,8 @@ "\n", "Scaling up a real cluster to a larger member population is then a CSV-level edit: add a row per\n", "additional member, fill in its centre and luminosity, save. The modeling script picks up the new rows\n", - "automatically and the scaling tier's free-parameter count stays at 1 (``b0_ref``, the reference\n", - "member's lens strength; the relation's exponents stay fixed at 0.5). Only luminosity *ratios* enter the\n", + "automatically and the scaling tier's free-parameter count stays at 1 (``sigma_ref``, the reference\n", + "member's fiducial velocity dispersion; the relation's exponents stay fixed). Only luminosity *ratios* enter the\n", "relation, so any consistent luminosity convention works \u2014 including converting from magnitudes.\n", "'''" ], diff --git a/notebooks/cluster/start_here.ipynb b/notebooks/cluster/start_here.ipynb index 6f8ba6b9c..34f2b41bd 100644 --- a/notebooks/cluster/start_here.ipynb +++ b/notebooks/cluster/start_here.ipynb @@ -327,16 +327,18 @@ "\n", " - **Main lens galaxies (2):** the two brightest core members (the BCG region galaxies),\n", " individually-modelled ``dPIEMassSph`` profiles with centre fixed to the observed light centres\n", - " and free ``ra``, ``rs``, ``b0``. **6 free parameters total.**\n", - "\n", - " - **Scaling-tier members (188):** ``dPIEMassSph`` profiles with centre fixed to the CSV centres. ``ra``,\n", - " ``rs`` and ``b0`` all derive from the reference-anchored relation used by Lenstool and standard in\n", - " published cluster analyses: ``ra = ra_ref * (L / L_ref) ** 0.5``, ``rs = rs_ref * (L / L_ref) ** 0.5``\n", - " and ``b0 = b0_ref * (L / L_ref) ** 0.5``, where ``L_ref`` is an explicit fixed reference luminosity\n", - " (Lenstool's ``mag0``), *not* the sample max. The exponent is fixed at the Faber-Jackson value\n", - " (b0 \u221d sigma\u00b2 and sigma \u221d L^(1/4) give b0 \u221d L^(1/2)) \u2014 only the normalization ``b0_ref``, the lens\n", - " strength of a reference-magnitude galaxy, is fitted. Our member luminosities are normalised to the\n", - " BCG's F160W flux, so ``L_ref = 1.0`` anchors the relation to the BCG itself.\n", + " and free ``sigma``, ``r_core``, ``r_cut`` \u2014 Lenstool's native dPIE parameters, so the posterior\n", + " reads like a Lenstool results table. **6 free parameters total.**\n", + "\n", + " - **Scaling-tier members (188):** ``dPIEMassSph`` profiles with centre fixed to the CSV centres.\n", + " ``sigma``, ``r_core`` and ``r_cut`` all derive from the reference-anchored relation used by Lenstool\n", + " and standard in published cluster analyses: ``sigma = sigma_ref * (L / L_ref) ** 0.25``,\n", + " ``r_core = r_core_ref * (L / L_ref) ** 0.5`` and ``r_cut = r_cut_ref * (L / L_ref) ** 0.5``, where\n", + " ``L_ref`` is an explicit fixed reference luminosity (Lenstool's ``mag0``), *not* the sample max. The\n", + " exponents are fixed at the Faber-Jackson values (L \u221d sigma^4 gives sigma \u221d L^(1/4); constant M/L gives\n", + " r_cut \u221d L^(1/2)) \u2014 only the normalization ``sigma_ref``, the fiducial velocity dispersion of a\n", + " reference-magnitude galaxy, is fitted. Our member luminosities are normalised to the BCG's F160W flux,\n", + " so ``L_ref = 1.0`` anchors the relation to the BCG itself.\n", " **1 free parameter total for the whole tier \u2014 independent of the number of members.**\n", "\n", " - **Host dark matter halo:** a standalone ``Galaxy`` carrying an ``NFWMCRLudlowSph`` halo with\n", @@ -365,7 +367,7 @@ "__Model__\n", "\n", "The model is composed below in four blocks: main-tier loop, host halo, source-tier loop, scaling-tier\n", - "loop (defining the shared ``b0_ref`` normalization once outside the loop). The four blocks are then\n", + "loop (defining the shared ``sigma_ref`` normalization once outside the loop). The four blocks are then\n", "bundled into a single ``af.Collection`` model that the analysis will receive." ] }, @@ -382,11 +384,16 @@ "\n", "galaxy_models = al.galaxy_af_models_from_csv_tables(mass_table, point_table)\n", "\n", - "# Main Lens Galaxies: free dPIE ra / rs / b0; centre stays fixed at the CSV value.\n", + "# Main Lens Galaxies: free dPIE sigma / r_core / r_cut; centre and redshifts stay\n", + "# fixed at the CSV values, and the cosmology constants H0 / Om0 are pinned (they are\n", + "# model *constants*, not parameters to sample \u2014 if left unset they would inherit the\n", + "# config's default priors and float).\n", "for name in (\"lens_0\", \"lens_1\"):\n", - " galaxy_models[name].mass.ra = af.UniformPrior(lower_limit=1.0, upper_limit=15.0)\n", - " galaxy_models[name].mass.rs = af.UniformPrior(lower_limit=5.0, upper_limit=40.0)\n", - " galaxy_models[name].mass.b0 = af.UniformPrior(lower_limit=0.1, upper_limit=10.0)\n", + " galaxy_models[name].mass.sigma = af.UniformPrior(lower_limit=50.0, upper_limit=600.0)\n", + " galaxy_models[name].mass.r_core = af.UniformPrior(lower_limit=1.0, upper_limit=15.0)\n", + " galaxy_models[name].mass.r_cut = af.UniformPrior(lower_limit=5.0, upper_limit=40.0)\n", + " galaxy_models[name].mass.H0 = 67.66\n", + " galaxy_models[name].mass.Om0 = 0.30966\n", "\n", "# Host Halo: free mass_at_200; centre + redshift_object + redshift_source fixed.\n", "galaxy_models[\"host_halo\"].dark.mass_at_200 = af.LogUniformPrior(\n", @@ -406,19 +413,21 @@ " mean=float(np.mean(positions[:, 1])), sigma=3.0\n", " )\n", "\n", - "# Scaling Tier (reference-anchored: b0_ref is the single shared free parameter, the\n", - "# lens strength of a galaxy at the reference magnitude; per-member ra, rs and b0\n", - "# derive from it with the exponent fixed at the Faber-Jackson value 0.5 \u2014 the\n", - "# Lenstool convention. The reference luminosity is an EXPLICIT FIXED constant\n", - "# (Lenstool's \"mag0\"); our member luminosities are normalised to the BCG's F160W\n", - "# flux, so L_ref = 1.0 anchors the relation to the BCG).\n", + "# Scaling Tier (reference-anchored: sigma_ref is the single shared free parameter,\n", + "# the fiducial velocity dispersion of a galaxy at the reference magnitude, in km/s;\n", + "# per-member sigma, r_core and r_cut derive from it with the exponents fixed at the\n", + "# Faber-Jackson values (sigma \u221d L^0.25, radii \u221d L^0.5) \u2014 the Lenstool potfile\n", + "# convention. The reference luminosity is an EXPLICIT FIXED constant (Lenstool's\n", + "# \"mag0\"); our member luminosities are normalised to the BCG's F160W flux, so\n", + "# L_ref = 1.0 anchors the relation to the BCG).\n", "\n", - "scaling_b0_ref = af.UniformPrior(lower_limit=0.0, upper_limit=1.0)\n", - "scaling_exponent = 0.5\n", + "scaling_sigma_ref = af.UniformPrior(lower_limit=0.0, upper_limit=300.0)\n", + "scaling_sigma_exponent = 0.25\n", + "scaling_radius_exponent = 0.5\n", "\n", "reference_luminosity = 1.0\n", - "scaling_ra_ref_fixed = 0.158\n", - "scaling_rs_ref_fixed = 15.8\n", + "scaling_r_core_ref_fixed = 0.158\n", + "scaling_r_cut_ref_fixed = 15.8\n", "\n", "scaling_galaxies_list = []\n", "for centre, luminosity in zip(\n", @@ -428,9 +437,13 @@ "\n", " mass = af.Model(al.mp.dPIEMassSph)\n", " mass.centre = tuple(centre)\n", - " mass.ra = scaling_ra_ref_fixed * luminosity_ratio**scaling_exponent\n", - " mass.rs = scaling_rs_ref_fixed * luminosity_ratio**scaling_exponent\n", - " mass.b0 = scaling_b0_ref * luminosity_ratio**scaling_exponent\n", + " mass.sigma = scaling_sigma_ref * luminosity_ratio**scaling_sigma_exponent\n", + " mass.r_core = scaling_r_core_ref_fixed * luminosity_ratio**scaling_radius_exponent\n", + " mass.r_cut = scaling_r_cut_ref_fixed * luminosity_ratio**scaling_radius_exponent\n", + " mass.redshift_object = redshift_lens\n", + " mass.redshift_source = max(source_redshifts)\n", + " mass.H0 = 67.66\n", + " mass.Om0 = 0.30966\n", "\n", " scaling_galaxies_list.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass))\n", "\n", diff --git a/notebooks/group/features/linear_light_profiles/slam.ipynb b/notebooks/group/features/linear_light_profiles/slam.ipynb index 9570c74fb..b2e35e84c 100644 --- a/notebooks/group/features/linear_light_profiles/slam.ipynb +++ b/notebooks/group/features/linear_light_profiles/slam.ipynb @@ -324,9 +324,21 @@ "\n", " extra_galaxies = af.Collection(extra_mass_models) if extra_mass_models else None\n", "\n", - " # --- scaling lens galaxy models (light fixed, shared luminosity scaling relation) ---\n", - " scaling_factor = af.UniformPrior(lower_limit=0.0, upper_limit=0.5)\n", - " scaling_relation = af.UniformPrior(lower_limit=0.0, upper_limit=2.0)\n", + " # --- scaling lens galaxy models (light fixed, reference-anchored scaling relation) ---\n", + " # Lenstool convention: einstein_radius = einstein_radius_ref * (L / L_ref)^0.5, exponent\n", + " # FIXED at the Faber-Jackson value; the single free parameter einstein_radius_ref is the\n", + " # Einstein radius of a galaxy as bright as the reference. The reference luminosity is the\n", + " # brightest MAIN lens galaxy (the BGG) \u2014 only its luminosity anchors the relation, its\n", + " # (free) mass is NOT coupled to the tier. See scripts/group/slam.py for the full rationale.\n", + " main_luminosity_list = [\n", + " abs(tracer.galaxies[i].bulge.luminosity_within_circle_from(radius=10.0))\n", + " / pixel_scale**2\n", + " for i in range(n_main)\n", + " ]\n", + " reference_luminosity = max(main_luminosity_list)\n", + "\n", + " einstein_radius_ref = af.UniformPrior(lower_limit=0.0, upper_limit=2.0)\n", + " scaling_exponent = 0.5\n", "\n", " scaling_mass_models = []\n", " for i in range(n_scaling):\n", @@ -340,7 +352,10 @@ " abs(galaxy_with_intensity.bulge.luminosity_within_circle_from(radius=10.0))\n", " / pixel_scale**2\n", " )\n", - " mass.einstein_radius = scaling_factor * total_luminosity**scaling_relation\n", + " mass.einstein_radius = (\n", + " einstein_radius_ref\n", + " * (total_luminosity / reference_luminosity) ** scaling_exponent\n", + " )\n", "\n", " scaling_mass_models.append(\n", " af.Model(\n", @@ -827,9 +842,18 @@ "\n", " extra_galaxies = af.Collection(extra_mass_models) if extra_mass_models else None\n", "\n", - " # --- scaling galaxies: fixed light, free shared scaling relation ---\n", - " scaling_factor = af.UniformPrior(lower_limit=0.0, upper_limit=0.5)\n", - " scaling_relation = af.UniformPrior(lower_limit=0.0, upper_limit=2.0)\n", + " # --- scaling galaxies: fixed light, reference-anchored scaling relation ---\n", + " # Same Lenstool convention as source_lp_1: free einstein_radius_ref, exponent fixed\n", + " # at 0.5, reference luminosity = the brightest main lens (BGG).\n", + " main_luminosity_list = [\n", + " abs(tracer.galaxies[i].bulge.luminosity_within_circle_from(radius=10.0))\n", + " / pixel_scale**2\n", + " for i in range(n_lenses)\n", + " ]\n", + " reference_luminosity = max(main_luminosity_list)\n", + "\n", + " einstein_radius_ref = af.UniformPrior(lower_limit=0.0, upper_limit=2.0)\n", + " scaling_exponent = 0.5\n", "\n", " scaling_mass_models = []\n", " for i in range(n_scaling):\n", @@ -843,7 +867,10 @@ " abs(galaxy_with_intensity.bulge.luminosity_within_circle_from(radius=10.0))\n", " / pixel_scale**2\n", " )\n", - " mass.einstein_radius = scaling_factor * total_luminosity**scaling_relation\n", + " mass.einstein_radius = (\n", + " einstein_radius_ref\n", + " * (total_luminosity / reference_luminosity) ** scaling_exponent\n", + " )\n", "\n", " scaling_mass_models.append(\n", " af.Model(\n", diff --git a/notebooks/group/features/scaling_relation/modeling.ipynb b/notebooks/group/features/scaling_relation/modeling.ipynb index cee97224c..b8006d4e0 100644 --- a/notebooks/group/features/scaling_relation/modeling.ipynb +++ b/notebooks/group/features/scaling_relation/modeling.ipynb @@ -385,9 +385,10 @@ "normalization-slope degeneracy. Only luminosity ratios enter, so the luminosity units are irrelevant; magnitude\n", "catalogues convert via ``L / L_ref = 10 ** (0.4 * (m_ref - m))``.\n", "\n", - "The dPIE-profile cluster-scale analogue \u2014 which also scales the core and truncation radii (``ra, rs \u221d L^0.5``,\n", - "mirroring Lenstool's r_core / r_cut scaling) \u2014 is ``scripts/cluster/modeling.py``. To free the exponent as a\n", - "systematics test, replace the fixed value with e.g. ``af.UniformPrior(lower_limit=0.0, upper_limit=1.0)``.\n", + "The dPIE-profile cluster-scale analogue \u2014 expressed directly in Lenstool's native parameters\n", + "(``sigma \u221d L^0.25``, ``r_core, r_cut \u221d L^0.5``, matching Lenstool's potfile scalings; equivalent to\n", + "einstein_radius \u221d L^0.5 since the lens strength goes as sigma\u00b2) \u2014 is ``scripts/cluster/modeling.py``. To free the\n", + "exponent as a systematics test, replace the fixed value with e.g. ``af.UniformPrior(lower_limit=0.0, upper_limit=1.0)``.\n", "\n", "Adding more scaling galaxies (e.g. by lengthening the centres + luminosity lists) does not add any free parameters\n", "to the model." diff --git a/notebooks/group/slam.ipynb b/notebooks/group/slam.ipynb index 3af681224..c2f8109e1 100644 --- a/notebooks/group/slam.ipynb +++ b/notebooks/group/slam.ipynb @@ -59,12 +59,18 @@ "\n", "Scaling galaxies are a larger ensemble of companions whose masses are constrained through a shared\n", "luminosity-to-mass scaling relation rather than being individually free. They each carry a free MGE\n", - "light profile, but their Einstein radii follow:\n", - "\n", - " einstein_radius = scaling_factor * total_luminosity^scaling_relation\n", - "\n", - "where `scaling_factor` and `scaling_relation` are two shared free parameters whose priors are\n", - "`UniformPrior(0, 0.5)` and `UniformPrior(0, 2)` respectively. This reduces the number of mass\n", + "light profile, but their Einstein radii follow the reference-anchored convention used by Lenstool\n", + "and essentially every published group- and cluster-scale analysis (Limousin et al. 2005;\n", + "Eliasdottir et al. 2007; Bergamini et al. 2019):\n", + "\n", + " einstein_radius = einstein_radius_ref * (luminosity / reference_luminosity)^0.5\n", + "\n", + "The exponent is FIXED at the Faber-Jackson value (einstein_radius \u221d sigma\u00b2 and sigma \u221d L^(1/4) give\n", + "einstein_radius \u221d L^(1/2)), avoiding the normalization-slope degeneracy of fitting both. The tier's\n", + "single shared free parameter is `einstein_radius_ref` \u2014 the Einstein radius of a galaxy as bright as\n", + "the reference \u2014 with the reference luminosity anchored to the brightest main lens galaxy (the BGG),\n", + "whose luminosity is fixed by the preceding light fit. Only the BGG's *luminosity* anchors the\n", + "relation; its (free) mass is deliberately not coupled to the tier. This reduces the number of mass\n", "parameters considerably when many companion galaxies are present. Scaling-galaxy models are stored\n", "in `model.scaling_galaxies`.\n", "\n", @@ -274,9 +280,10 @@ "\n", "Multiple main-lens galaxies each get an `Isothermal` mass; only `lens_0` carries an\n", "`ExternalShear`. Extra-galaxy Einstein radii are bounded by a luminosity-derived prior\n", - "(`min(5 * 0.5 * L^0.6, 5.0)`). Scaling galaxies share two free parameters,\n", - "`scaling_factor` and `scaling_relation`, so their masses follow\n", - "`einstein_radius = scaling_factor * luminosity^scaling_relation`." + "(`min(5 * 0.5 * L^0.6, 5.0)`). Scaling galaxies share one free parameter,\n", + "`einstein_radius_ref`, so their masses follow the reference-anchored Faber-Jackson relation\n", + "`einstein_radius = einstein_radius_ref * (luminosity / reference_luminosity)^0.5`, with the\n", + "reference luminosity fixed to the brightest main lens (BGG)." ] }, { @@ -377,9 +384,30 @@ "\n", " extra_galaxies = af.Collection(extra_mass_models) if extra_mass_models else None\n", "\n", - " # --- scaling lens galaxy models (light fixed, shared luminosity scaling relation) ---\n", - " scaling_factor = af.UniformPrior(lower_limit=0.0, upper_limit=0.5)\n", - " scaling_relation = af.UniformPrior(lower_limit=0.0, upper_limit=2.0)\n", + " # --- scaling lens galaxy models (light fixed, reference-anchored scaling relation) ---\n", + " # Lenstool convention: einstein_radius = einstein_radius_ref * (L / L_ref)^0.5 with the\n", + " # exponent FIXED at the Faber-Jackson value (einstein_radius \u221d sigma\u00b2 and sigma \u221d L^(1/4)\n", + " # give einstein_radius \u221d L^(1/2)). The tier's single free parameter is\n", + " # einstein_radius_ref \u2014 the Einstein radius of a galaxy as bright as the reference \u2014 a\n", + " # physically interpretable normalization with an easy-to-motivate prior range, avoiding\n", + " # the normalization-slope degeneracy of fitting both. The reference luminosity is the\n", + " # brightest MAIN lens galaxy (the BGG) \u2014 an explicit anchor fixed by the preceding light\n", + " # fit and invariant to which galaxies are in the tier. Only the BGG's *luminosity* is the\n", + " # anchor; its (free) mass is deliberately NOT coupled to the tier.\n", + " main_luminosity_list = [\n", + " np.sum(\n", + " [\n", + " 2 * np.pi * g.sigma**2 / g.axis_ratio() * g.intensity\n", + " for g in tracer.galaxies[i].bulge.profile_list\n", + " ]\n", + " )\n", + " / pixel_scale**2\n", + " for i in range(n_main)\n", + " ]\n", + " reference_luminosity = max(main_luminosity_list)\n", + "\n", + " einstein_radius_ref = af.UniformPrior(lower_limit=0.0, upper_limit=2.0)\n", + " scaling_exponent = 0.5\n", "\n", " scaling_mass_models = []\n", " for i in range(n_scaling):\n", @@ -394,7 +422,10 @@ " for g in tracer.galaxies[n_main + n_extra + i].bulge.profile_list\n", " ]\n", " total_luminosity = np.sum(luminosity_per_gaussian_list) / pixel_scale**2\n", - " mass.einstein_radius = scaling_factor * total_luminosity**scaling_relation\n", + " mass.einstein_radius = (\n", + " einstein_radius_ref\n", + " * (total_luminosity / reference_luminosity) ** scaling_exponent\n", + " )\n", "\n", " scaling_mass_models.append(\n", " af.Model(\n", @@ -894,9 +925,25 @@ "\n", " extra_galaxies = af.Collection(extra_mass_models) if extra_mass_models else None\n", "\n", - " # --- scaling galaxies: fixed light, free shared scaling relation ---\n", - " scaling_factor = af.UniformPrior(lower_limit=0.0, upper_limit=0.5)\n", - " scaling_relation = af.UniformPrior(lower_limit=0.0, upper_limit=2.0)\n", + " # --- scaling galaxies: fixed light, reference-anchored scaling relation ---\n", + " # Same Lenstool convention as source_lp_1: free einstein_radius_ref, exponent fixed\n", + " # at 0.5, reference luminosity = the brightest main lens (BGG) from the light[1] fit.\n", + " main_luminosity_list = [\n", + " np.sum(\n", + " [\n", + " 2 * np.pi * g.sigma**2 / g.axis_ratio() * g.intensity\n", + " for g in getattr(\n", + " light_result.instance.galaxies, f\"lens_{i}\"\n", + " ).bulge.profile_list\n", + " ]\n", + " )\n", + " / pixel_scale**2\n", + " for i in range(n_lenses)\n", + " ]\n", + " reference_luminosity = max(main_luminosity_list)\n", + "\n", + " einstein_radius_ref = af.UniformPrior(lower_limit=0.0, upper_limit=2.0)\n", + " scaling_exponent = 0.5\n", "\n", " scaling_mass_models = []\n", " for i in range(n_scaling):\n", @@ -911,7 +958,10 @@ " for g in light_scaling.bulge.profile_list\n", " ]\n", " total_luminosity = np.sum(luminosity_per_gaussian_list) / pixel_scale**2\n", - " mass.einstein_radius = scaling_factor * total_luminosity**scaling_relation\n", + " mass.einstein_radius = (\n", + " einstein_radius_ref\n", + " * (total_luminosity / reference_luminosity) ** scaling_exponent\n", + " )\n", "\n", " scaling_mass_models.append(\n", " af.Model(\n", diff --git a/notebooks/guides/advanced/potential_correction.ipynb b/notebooks/guides/advanced/potential_correction.ipynb new file mode 100644 index 000000000..06a526b87 --- /dev/null +++ b/notebooks/guides/advanced/potential_correction.ipynb @@ -0,0 +1,481 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Potential Corrections (Gravitational Imaging)\n", + "=============================================\n", + "\n", + "Parametric lens mass models (e.g. the `PowerLaw` used throughout the workspace) are smooth: their convergence and\n", + "deflection angles vary gradually across the image plane. Real lens galaxies are not perfectly smooth. They contain\n", + "substructure \u2014 dark matter subhaloes, globular clusters, satellite galaxies \u2014 and their large-scale mass\n", + "distributions can depart from ellipsoidal symmetry in ways no simple parametric profile captures.\n", + "\n", + "Gravitational imaging, also called the potential correction method, reconstructs these departures directly. Small\n", + "pixelized corrections $\\delta\\psi$ to the lensing potential are defined on a coarse regular mesh over the image\n", + "plane, and solved for linearly alongside the pixelized source, by maximising the Bayesian evidence. Wherever the\n", + "smooth model fails to fit the lensed arcs, the corrections absorb the discrepancy \u2014 and their curvature\n", + "$\\delta\\kappa = \\frac{1}{2}\\nabla^2 \\delta\\psi$ maps the missing (or excess) convergence, revealing substructure\n", + "without assuming its form.\n", + "\n", + "This technique detected a $\\sim 10^6 \\, M_\\odot$ object at cosmological distance in the VLBI imaging of\n", + "JVAS B1938+666 (Powell et al. 2025, Nature Astronomy 9, 1714; Vegetti et al. 2026, Nature Astronomy), the smallest\n", + "dark object ever found by lensing. The implementation in **PyAutoLens** (`al.pc`) is ported from the\n", + "`potential_correction` package of Cao et al. 2025 (https://github.com/caoxiaoyue/lensing_potential_correction).\n", + "If you use it in your research, please cite Cao et al. 2025 \u2014 citation materials are at\n", + "https://github.com/caoxiaoyue/potential_correction_paper.\n", + "\n", + "__Overview__\n", + "\n", + "The method models the observed image $d$ as the smooth-model image plus a linear response to the corrections:\n", + "\n", + "$\\delta d = - B \\, D_s \\, D_\\psi \\, \\delta\\psi$\n", + "\n", + "where $B$ is the PSF blur matrix, $D_s$ holds the source's brightness gradients at the ray-traced position of every\n", + "image pixel, and $D_\\psi$ interpolates the coarse $\\delta\\psi$ mesh onto the image grid and takes its spatial\n", + "gradients. Intuitively: a small change to the potential deflects a ray slightly, which samples the source at a\n", + "slightly different position, which changes the observed brightness in proportion to the source's local gradient.\n", + "\n", + "Both the source and $\\delta\\psi$ are regularized, and the strengths of both regularizations are set objectively by\n", + "maximising the Bayesian evidence \u2014 no manual tuning, following Cao et al. 2025 (building on Koopmans 2005;\n", + "Vegetti & Koopmans 2009; Suyu et al. 2009; Vernardos & Koopmans 2022).\n", + "\n", + "__Contents__\n", + "\n", + "- **Simulate:** Simulate strong lens imaging whose mass distribution contains a dark subhalo the smooth model omits.\n", + "- **Smooth Model Fit:** Fit the data with the smooth (subhalo-free) mass model, whose residuals localize the subhalo.\n", + "- **Joint Fit:** Reconstruct the source and potential corrections jointly with `FitDpsiSrcImaging`.\n", + "- **Dkappa Map:** Convert the corrections to a convergence map and compare to the true subhalo.\n", + "- **Iterative Fit:** Refine the corrections with the iterative Levenberg-Marquardt engine `IterFitDpsiSrcImaging`.\n", + "- **Evidence Sampling:** Sample the regularization hyper-parameters with a non-linear search via `DpsiSrcInvAnalysis`." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Google Colab Setup__\n", + "\n", + "This cell sets up the environment when the notebook is run on Google Colab: it installs the\n", + "required PyAuto packages, clones the workspace (configuration files and example datasets) and\n", + "points the configuration at it. If you are running the notebook elsewhere (e.g. locally via\n", + "your own installation) it does nothing, and you can run it safely.\n", + "\n", + "Colab tip: model-fits run much faster on a GPU \u2014 enable one via \"Runtime\" -> \"Change runtime\n", + "type\" -> \"Hardware accelerator\" before running the notebook." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "try:\n", + " import google.colab\n", + " import subprocess\n", + " import sys\n", + "\n", + " subprocess.check_call(\n", + " [sys.executable, \"-m\", \"pip\", \"install\", \"autoconf\", \"--no-deps\"]\n", + " )\n", + "except ImportError:\n", + " pass\n", + "\n", + "from autoconf import setup_colab\n", + "\n", + "setup_colab.setup(\"autolens\")" + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "\n", + "# %matplotlib inline\n", + "# from pyprojroot import here\n", + "# workspace_path = str(here())\n", + "# %cd $workspace_path\n", + "# print(f\"Working Directory has been set to `{workspace_path}`\")\n", + "\n", + "import numpy as np\n", + "from os import path\n", + "\n", + "import autofit as af\n", + "import autolens as al\n", + "import autolens.plot as aplt" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Simulate__\n", + "\n", + "We simulate imaging of a strong lens whose mass model is an `Isothermal` plus a $10^{10} \\, M_\\odot$ NFW dark\n", + "subhalo sitting right on the Einstein ring, where its imprint on the lensed arcs is strongest.\n", + "\n", + "The source is a compact double-knot galaxy, giving the arcs the sharp brightness gradients that gravitational\n", + "imaging leverages (the response to $\\delta\\psi$ scales with the source gradient)." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "grid = al.Grid2D.uniform(shape_native=(120, 120), pixel_scales=0.05, over_sample_size=4)\n", + "\n", + "psf = al.Convolver.from_gaussian(shape_native=(11, 11), sigma=0.05, pixel_scales=0.05)\n", + "\n", + "simulator = al.SimulatorImaging(\n", + " exposure_time=840.0,\n", + " psf=psf,\n", + " background_sky_level=0.1,\n", + " add_poisson_noise_to_data=True,\n", + " noise_seed=1,\n", + ")\n", + "\n", + "lens_true = al.Galaxy(\n", + " redshift=0.2,\n", + " mass=al.mp.Isothermal(\n", + " centre=(0.0, 0.0),\n", + " einstein_radius=1.4,\n", + " ell_comps=al.convert.ell_comps_from(axis_ratio=0.9, angle=0.0),\n", + " ),\n", + " subhalo=al.mp.NFWMCRLudlowSph(\n", + " centre=(1.41, 0.0),\n", + " mass_at_200=1.0e10,\n", + " redshift_object=0.2,\n", + " redshift_source=0.6,\n", + " ),\n", + ")\n", + "\n", + "source_true = al.Galaxy(\n", + " redshift=0.6,\n", + " bulge0=al.lp.Gaussian(\n", + " centre=(0.0, 0.0),\n", + " ell_comps=al.convert.ell_comps_from(axis_ratio=0.6, angle=45.0),\n", + " intensity=5.0,\n", + " sigma=0.15,\n", + " ),\n", + " bulge1=al.lp.Gaussian(\n", + " centre=(0.0, 0.4),\n", + " ell_comps=al.convert.ell_comps_from(axis_ratio=0.4, angle=135.0),\n", + " intensity=3.0,\n", + " sigma=0.1,\n", + " ),\n", + ")\n", + "\n", + "tracer_true = al.Tracer(galaxies=[lens_true, source_true])\n", + "dataset = simulator.via_tracer_from(tracer=tracer_true, grid=grid)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Arc Mask__\n", + "\n", + "Gravitational imaging only constrains the potential where the lensed arcs are: the corrections respond to the data\n", + "through the source gradients, which vanish off the arcs. We therefore mask the data to an arc-tracing region using\n", + "`al.pc.util.arc_mask_from`, which thresholds the signal-to-noise map, drops small disconnected islands, dilates the\n", + "result and cleans it so every unmasked pixel supports the finite-difference derivative operators." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "mask_array = al.pc.util.arc_mask_from(\n", + " np.asarray(dataset.signal_to_noise_map.native), threshold=3.0, ignore_size=25, ext_size=5\n", + ")\n", + "mask = al.Mask2D(mask=mask_array, pixel_scales=dataset.pixel_scales)\n", + "masked_imaging = dataset.apply_mask(mask=mask)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Smooth Model Fit__\n", + "\n", + "The starting point of a potential-correction analysis is the best smooth-model fit \u2014 in a real analysis, the result\n", + "of a standard lens-modeling pipeline. Its residuals concentrate near the subhalo, because the smooth model cannot\n", + "bend the arcs the way the subhalo does.\n", + "\n", + "Here we cheat and use the true smooth mass model (the `Isothermal` without the subhalo) and the true source, which\n", + "isolates the technique itself; a real analysis would use the maximum-likelihood parametric model.\n", + "\n", + "The source enters the correction operator through an `al.pc.SrcFactory`, which evaluates the source's brightness and\n", + "gradients at arbitrary source-plane positions. For a parametric source we use `AnalyticSrcFactory`; a pixelized\n", + "reconstruction from a previous fit would use `PixSrcFactoryITP`." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "lens_smooth = al.Galaxy(redshift=0.2, mass=lens_true.mass)\n", + "source_start = al.pc.AnalyticSrcFactory(source_galaxy=source_true)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Joint Fit__\n", + "\n", + "`FitDpsiSrcImaging` inverts the image jointly for the pixelized source and the potential corrections, fully\n", + "accounting for their covariance. Its components are:\n", + "\n", + "- `DpsiPixelization`: the corrections' mesh (`RegularDpsiMesh(factor=2)` \u2014 a mesh twice as coarse as the data grid)\n", + " and their regularization. `aa.reg.MaternKernel` (with `nu=2.5`) flexibly permits both localised and extended\n", + " perturbations (the key advance of Cao et al. 2025 over curvature-only schemes); `aa.reg.CurvatureMask` and\n", + " `aa.reg.FourthOrderMask` are the classic alternatives.\n", + "\n", + "- The source pixelization: any standard PyAutoLens pixelization; here a k-nearest-neighbour mesh distributed by an\n", + " `Overlay` image mesh." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "dpsi_pixelization = al.pc.DpsiPixelization(\n", + " mesh=al.pc.RegularDpsiMesh(factor=2),\n", + " regularization=al.reg.MaternKernel(coefficient=2000.0, scale=4.0, nu=2.5),\n", + ")\n", + "\n", + "grid_slim = masked_imaging.grid.slim\n", + "source_shape = (\n", + " int(float(grid_slim[:, 0].max() - grid_slim[:, 0].min()) / 0.05 / 2.0),\n", + " int(float(grid_slim[:, 1].max() - grid_slim[:, 1].min()) / 0.05 / 2.0),\n", + ")\n", + "src_pixelization = al.Pixelization(\n", + " mesh=al.mesh.KNearestNeighbor(pixels=int(np.prod(source_shape))),\n", + " regularization=al.reg.Constant(coefficient=3.8),\n", + ")\n", + "src_image_mesh = al.image_mesh.Overlay(shape=source_shape)\n", + "\n", + "fit = al.pc.FitDpsiSrcImaging(\n", + " masked_imaging=masked_imaging,\n", + " lens_start=lens_smooth,\n", + " source_start=source_start,\n", + " dpsi_pixelization=dpsi_pixelization,\n", + " src_pixelization=src_pixelization,\n", + " src_image_mesh=src_image_mesh,\n", + ")\n", + "\n", + "print(f\"joint source + dpsi log evidence = {fit.log_evidence:.4e}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Dkappa Map__\n", + "\n", + "The reconstructed corrections live on the coarse dpsi mesh (`fit.best_fit_dpsi`). Their physical meaning is clearest\n", + "as a convergence correction, obtained by applying the mesh's Laplacian operator:\n", + "\n", + "$\\delta\\kappa = \\frac{1}{2} \\nabla^2 \\delta\\psi$\n", + "\n", + "A dark subhalo the smooth model omits appears as a positive $\\delta\\kappa$ peak at its position \u2014 compare the map\n", + "below to the true subhalo at (y, x) = (1.41\", 0.0\").\n", + "\n", + "The nine-panel summary figure shows the data, model, residuals, the $\\delta\\psi$ and $\\delta\\kappa$ maps and the\n", + "source reconstruction." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "dkappa = fit.pair_dpsi_data_obj.hamiltonian_dpsi @ fit.best_fit_dpsi\n", + "\n", + "peak = np.argmax(dkappa)\n", + "print(\n", + " f\"dkappa peak at (y, x) = ({fit.pair_dpsi_data_obj.ygrid_dpsi_1d[peak]:.2f}, \"\n", + " f\"{fit.pair_dpsi_data_obj.xgrid_dpsi_1d[peak]:.2f}) \u2014 true subhalo at (1.41, 0.00)\"\n", + ")\n", + "\n", + "al.pc.visualize.show_fit_dpsi_src(fit, output=path.join(\"output\", \"potential_correction_joint_fit.png\"))" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Iterative Fit__\n", + "\n", + "The joint fit above linearizes the corrections around the smooth model once. The iterative engine\n", + "`IterFitDpsiSrcImaging` goes further, following the approach of Powell et al. 2025: it optimizes the combined\n", + "state [source | dpsi] with a Levenberg-Marquardt loop, and after every accepted step **re-ray-traces the image grid\n", + "through the corrected lens** (the smooth model plus an `InputPotential` mass profile built from the current\n", + "corrections). The corrections thus feed back into the source mapping, capturing compact perturbers more faithfully\n", + "than a single linearization.\n", + "\n", + "`gauge_constraints=True` removes the degeneracies inherent to potential corrections \u2014 a constant $\\delta\\psi$ and\n", + "linear gradients (which just shift the source) are unconstrained by the data \u2014 by enforcing\n", + "$\\langle\\delta\\psi, 1\\rangle = \\langle\\delta\\psi, x\\rangle = \\langle\\delta\\psi, y\\rangle = 0$ at every step.\n", + "\n", + "The dense linear algebra runs through `al.pc.dense_util`, whose kernels follow the PyAuto `xp` convention: pass\n", + "`xp=jax.numpy` to `solve_joint_optimization` to run them under JAX on an accelerator; the default is numpy." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "iter_fit = al.pc.IterFitDpsiSrcImaging(\n", + " masked_imaging=masked_imaging,\n", + " lens_start=lens_smooth,\n", + " dpsi_pixelization=dpsi_pixelization,\n", + " src_pixelization=src_pixelization,\n", + " src_image_mesh=src_image_mesh,\n", + " gauge_constraints=True,\n", + " n_iter=5,\n", + ")\n", + "\n", + "s_opt, dpsi_opt = iter_fit.solve_joint_optimization()\n", + "print(f\"iterative Laplace log evidence = {iter_fit.log_evidence():.4e}\")\n", + "\n", + "dkappa_iter = iter_fit.pair_dpsi_data_obj.hamiltonian_dpsi @ dpsi_opt\n", + "peak = np.argmax(dkappa_iter)\n", + "print(\n", + " f\"iterative dkappa peak at (y, x) = ({iter_fit.pair_dpsi_data_obj.ygrid_dpsi_1d[peak]:.2f}, \"\n", + " f\"{iter_fit.pair_dpsi_data_obj.xgrid_dpsi_1d[peak]:.2f})\"\n", + ")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Evidence Sampling__\n", + "\n", + "In a real analysis the regularization hyper-parameters (the coefficients and scales above) are not known. They are\n", + "sampled with a non-linear search, using the inversion's Bayesian evidence as the likelihood, via the\n", + "`DpsiSrcInvAnalysis` (one-shot) or `IterDpsiSrcInvAnalysis` (iterative) analysis classes:\n", + "\n", + " dpsi_model = af.Model(\n", + " al.pc.DpsiSrcPixelization,\n", + " dpsi_pixelization=af.Model(\n", + " al.pc.DpsiPixelization,\n", + " mesh=al.pc.RegularDpsiMesh(factor=2),\n", + " regularization=af.Model(al.reg.MaternKernel, nu=2.5),\n", + " ),\n", + " src_pixelization=src_pixelization,\n", + " )\n", + "\n", + " analysis = al.pc.DpsiSrcInvAnalysis(\n", + " masked_imaging=masked_imaging,\n", + " lens_start=lens_smooth,\n", + " source_start=source_start,\n", + " src_image_mesh=src_image_mesh,\n", + " )\n", + "\n", + " search = af.Nautilus(name=\"potential_correction\", n_live=100)\n", + " result = search.fit(model=dpsi_model, analysis=analysis)\n", + "\n", + "This recovers both localised subhaloes and extended perturbations (e.g. Gaussian-random-field departures \u2014 see\n", + "`ag.mp.GaussianRandomField` for simulating them) with the regularization set objectively by the data.\n", + "\n", + "__Interferometer__\n", + "\n", + "Gravitational imaging's benchmark detections are radio/VLBI (the B1938+666 results above), and `al.pc` supports\n", + "`Interferometer` datasets in visibility space through the **sparse-operator (w-tilde) route**, whose cost scales\n", + "with real-space mask pixels independent of the visibility count. The certified recipe (validated on realistic\n", + "earth-rotation-synthesis coverage at ~10^4 visibilities, where both engines recover a simulated\n", + "$10^{10} M_\\odot$ subhalo at ~9$\\sigma$ with the dkappa peak 0.34\" from the truth) is:\n", + "\n", + " # the real-space mask stays a filled circle (it defines the FFT extent);\n", + " # the corrections are restricted to an arc-tracing sub-mask\n", + " dataset = dataset.apply_sparse_operator() # precompute the w-tilde operator (cache it to disk!)\n", + "\n", + " fit = al.pc.FitDpsiSrcInterferometer(\n", + " dataset=dataset,\n", + " lens_start=lens_smooth,\n", + " source_start=source_start,\n", + " dpsi_pixelization=dpsi_pixelization,\n", + " src_pixelization=src_pixelization,\n", + " dpsi_mask=arc_dpsi_mask,\n", + " )\n", + " evidence = fit.log_evidence # one-shot joint inversion (sparse route)\n", + "\n", + " iter_fit = al.pc.IterFitDpsiSrcInterferometer(\n", + " dataset=dataset,\n", + " lens_start=lens_smooth,\n", + " dpsi_pixelization=dpsi_pixelization,\n", + " src_pixelization=src_pixelization,\n", + " dpsi_mask=arc_dpsi_mask,\n", + " gauge_constraints=True,\n", + " reg_optimize_every=1, # evidence-control the regularizations each step\n", + " )\n", + " s_opt, dpsi_opt = iter_fit.solve_joint_optimization(\n", + " x0=fit.src_dpsi_slim, # warm-start from the one-shot solution\n", + " )\n", + "\n", + "Two practices matter, both mirroring the published pipelines: **warm-start the iterative engine from the\n", + "one-shot solution** (it then refines inside the right basin instead of searching from zero), and let the\n", + "**evidence control the regularization strengths** (`reg_optimize_every`) rather than fixing them \u2014 the Bayesian\n", + "evidence consistently ranks over-fit or degenerate solutions below the recovering ones, so evidence-driven\n", + "selection self-protects. Always check first that your smooth-model source fit reaches chi-squared per degree of\n", + "freedom near ~1-3: outside that regime the corrections absorb source-model error rather than real structure.\n", + "\n", + "__Wrap Up__\n", + "\n", + "Gravitational imaging turns the residuals of a smooth lens model into a map of the missing mass. The `al.pc`\n", + "subpackage provides the one-shot joint inversion, the iterative LM engine, and the analysis classes to sample their\n", + "hyper-parameters \u2014 with all linear algebra available under numpy or JAX through the `xp` API.\n", + "\n", + "If you use this functionality, please cite Cao et al. 2025\n", + "(https://github.com/caoxiaoyue/potential_correction_paper) alongside PyAutoLens." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [], + "outputs": [], + "execution_count": null + } + ], + "metadata": { + "anaconda-cloud": {}, + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.1" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} \ No newline at end of file diff --git a/notebooks/guides/profiles/mass.ipynb b/notebooks/guides/profiles/mass.ipynb index 8d5665395..45ed08f2d 100644 --- a/notebooks/guides/profiles/mass.ipynb +++ b/notebooks/guides/profiles/mass.ipynb @@ -191,11 +191,16 @@ "power_law_broken = al.mp.PowerLawBroken()\n", "power_law_broken_sph = al.mp.PowerLawBrokenSph()\n", "\n", - "# Pseudo-isothermal family (mass and potential parameterisations)\n", - "# Note: dPIEMass with default ell_comps=(0,0) triggers a divide-by-zero in the complex-plane\n", - "# formula; we use a small ellipticity here so the survey constructions succeed cleanly.\n", - "dpie_mass = al.mp.dPIEMass(ell_comps=(0.05, 0.0))\n", + "# Pseudo-isothermal family. The default dPIEMass / dPIEMassSph are parameterized in\n", + "# Lenstool's native convention (ellipticity, angle_pos, sigma = fiducial v_disp in km/s,\n", + "# r_core, r_cut, plus the redshifts entering the D_LS/D_S normalization) \u2014 the same numbers\n", + "# that appear in published cluster papers' results tables. The internal (ra, rs, b0)\n", + "# parameterization is the non-standard dPIEMassB0 / dPIEMassB0Sph.\n", + "dpie_mass = al.mp.dPIEMass(ellipticity=0.1)\n", "dpie_mass_sph = al.mp.dPIEMassSph()\n", + "dpie_mass_b0 = al.mp.dPIEMassB0(ell_comps=(0.05, 0.0))\n", + "# Note: PIEMass with ell_comps=(0,0) triggers a divide-by-zero in the complex-plane\n", + "# formula; we use a small ellipticity here so the survey constructions succeed cleanly.\n", "pie_mass = al.mp.PIEMass(ell_comps=(0.05, 0.0))\n", "dpie_potential = al.mp.dPIEPotential()\n", "dpie_potential_sph = al.mp.dPIEPotentialSph()\n", @@ -784,9 +789,11 @@ "aplt.plot_array(\n", " array=al.mp.dPIEMass(\n", " centre=(0.0, 0.0),\n", - " ell_comps=al.convert.ell_comps_from(axis_ratio=0.9, angle=45.0),\n", - " ra=0.1,\n", - " b0=0.5,\n", + " ellipticity=0.1,\n", + " angle_pos=45.0,\n", + " sigma=200.0,\n", + " r_core=0.1,\n", + " r_cut=20.0,\n", " ).convergence_2d_from(grid=grid),\n", " title=\"dPIEMass Convergence\",\n", ")\n", diff --git a/notebooks/imaging/features/potential_correction/likelihood_function.ipynb b/notebooks/imaging/features/potential_correction/likelihood_function.ipynb new file mode 100644 index 000000000..3a5e49c69 --- /dev/null +++ b/notebooks/imaging/features/potential_correction/likelihood_function.ipynb @@ -0,0 +1,643 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Log Likelihood Function: Potential Correction (Gravitational Imaging)__\n", + "\n", + "This script provides a step-by-step guide of the **PyAutoLens** potential-correction `log_likelihood_function`\n", + "(the Bayesian evidence of `al.pc.FitDpsiSrcImaging`), which jointly inverts `Imaging` data for a pixelized source\n", + "and pixelized corrections $\\delta\\psi$ to the lensing potential.\n", + "\n", + "This script has the following aims:\n", + "\n", + " - To provide a resource that authors can include in papers using **PyAutoLens**, so that readers can understand the\n", + " likelihood function (including references to the previous literature from which it is defined) without having to\n", + " write large quantities of text and equations.\n", + "\n", + " - To make gravitational imaging less of a \"black-box\" to users: every operator of the correction formalism\n", + " ($B$, $D_s$, $D_\\psi$) is built explicitly as a numpy array below.\n", + "\n", + "If you use the potential-correction functionality in your research, please cite Cao et al. 2025, from whose\n", + "`potential_correction` package (https://github.com/caoxiaoyue/lensing_potential_correction) the implementation is\n", + "ported; citation materials are provided at https://github.com/caoxiaoyue/potential_correction_paper. The method\n", + "builds on Koopmans 2005, Suyu et al. 2009, Vegetti & Koopmans 2009 and Vernardos & Koopmans 2022, and is the\n", + "image-plane analogue of the methodology behind the B1938+666 detections of Powell et al. 2025 (Nature Astronomy 9,\n", + "1714) and Vegetti et al. 2026.\n", + "\n", + "__Contents__\n", + "\n", + "- **Simplifications:** The choices made to keep this walkthrough small and explicit.\n", + "- **Prerequisites:** The likelihood functions this one builds on.\n", + "- **Dataset:** Simulate the imaging dataset fitted: a lens whose true mass contains a dark subhalo.\n", + "- **Arc Mask:** Mask the data to the arc region where the corrections are constrained.\n", + "- **Smooth Starting Model:** The smooth lens model and source the corrections perturb.\n", + "- **Dpsi Mesh:** The coarse rectangular mesh the corrections are defined on, paired to the data grid.\n", + "- **Dpsi Gradient Operator:** The sparse operator taking mesh corrections to their image-plane gradients.\n", + "- **Source Gradients:** The source's brightness gradients at the ray-traced position of every image pixel.\n", + "- **PSF Blur Matrix:** The explicit matrix form of the PSF convolution.\n", + "- **Dpsi Mapping Matrix:** Combining the three operators into the linear response -B D_s D_psi.\n", + "- **Source Inversion Blocks:** The standard pixelized source inversion supplying the source blocks.\n", + "- **Joint System:** The block mapping matrix and block-diagonal regularization of the joint inversion.\n", + "- **Solve:** Solving the joint linear system for the source and the corrections.\n", + "- **Dkappa Map:** Converting the solved corrections into a convergence-correction map.\n", + "- **Evidence Terms:** The five terms of the Bayesian evidence, computed explicitly.\n", + "- **Verification:** The same number from `al.pc.FitDpsiSrcImaging.log_evidence`.\n", + "- **JAX / xp:** The same evidence through the `al.pc.dense_util` xp-API kernels.\n", + "- **Wrap Up:** Summary and next steps.\n", + "\n", + "__Simplifications__\n", + "\n", + "To keep every array small enough to inspect, this example uses a low-resolution simulated dataset and a\n", + "`RectangularUniform` source mesh (uniform source pixels are simpler to reason about than adaptive meshes). The\n", + "corrections are regularized with the `MaternKernel` scheme (`nu=2.5`) used in science analyses \u2014 the key result of\n", + "Cao et al. 2025 is that the Matern family recovers localised perturbers where classic curvature-penalty schemes\n", + "(`aa.reg.CurvatureMask`, also usable here) smear them out. The smooth starting model and source are the truths\n", + "used in the simulation, isolating the correction formalism itself.\n", + "\n", + "__Prerequisites__\n", + "\n", + "- `imaging/likelihood_function.ipynb` \u2014 the standard imaging likelihood.\n", + "- `imaging/features/pixelization/likelihood_function.ipynb` \u2014 the pixelized-source evidence, whose terms reappear here.\n", + "- `guides/advanced/potential_correction.ipynb` \u2014 the user-facing overview of the technique." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Google Colab Setup__\n", + "\n", + "This cell sets up the environment when the notebook is run on Google Colab: it installs the\n", + "required PyAuto packages, clones the workspace (configuration files and example datasets) and\n", + "points the configuration at it. If you are running the notebook elsewhere (e.g. locally via\n", + "your own installation) it does nothing, and you can run it safely.\n", + "\n", + "Colab tip: model-fits run much faster on a GPU \u2014 enable one via \"Runtime\" -> \"Change runtime\n", + "type\" -> \"Hardware accelerator\" before running the notebook." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "try:\n", + " import google.colab\n", + " import subprocess\n", + " import sys\n", + "\n", + " subprocess.check_call(\n", + " [sys.executable, \"-m\", \"pip\", \"install\", \"autoconf\", \"--no-deps\"]\n", + " )\n", + "except ImportError:\n", + " pass\n", + "\n", + "from autoconf import setup_colab\n", + "\n", + "setup_colab.setup(\"autolens\")" + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "\n", + "from autoconf import jax_wrapper # Sets JAX environment before other imports\n", + "\n", + "from autoconf import setup_notebook; setup_notebook()\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from scipy.sparse import block_diag\n", + "\n", + "import autolens as al\n", + "import autolens.plot as aplt" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Dataset__\n", + "\n", + "We simulate the dataset in-memory (seeded, so this script is fully reproducible): an `IsothermalSph` lens whose\n", + "true mass also contains a $10^{10} M_\\odot$ NFW subhalo on the Einstein ring, lensing a compact Gaussian source.\n", + "The subhalo is what the corrections will recover." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "grid = al.Grid2D.uniform(shape_native=(80, 80), pixel_scales=0.08, over_sample_size=4)\n", + "psf = al.Convolver.from_gaussian(shape_native=(9, 9), sigma=0.08, pixel_scales=0.08)\n", + "\n", + "simulator = al.SimulatorImaging(\n", + " exposure_time=840.0,\n", + " psf=psf,\n", + " background_sky_level=0.1,\n", + " add_poisson_noise_to_data=True,\n", + " noise_seed=1,\n", + ")\n", + "\n", + "lens_true = al.Galaxy(\n", + " redshift=0.2,\n", + " mass=al.mp.IsothermalSph(centre=(0.0, 0.0), einstein_radius=1.4),\n", + " subhalo=al.mp.NFWMCRLudlowSph(\n", + " centre=(1.41, 0.0), mass_at_200=1.0e10, redshift_object=0.2, redshift_source=0.6\n", + " ),\n", + ")\n", + "source_true = al.Galaxy(\n", + " redshift=0.6,\n", + " bulge=al.lp.Gaussian(centre=(0.0, 0.0), intensity=5.0, sigma=0.2),\n", + ")\n", + "\n", + "dataset = simulator.via_tracer_from(\n", + " tracer=al.Tracer(galaxies=[lens_true, source_true]), grid=grid\n", + ")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Arc Mask__\n", + "\n", + "The corrections respond to the data only through the source's brightness gradients, which vanish away from the\n", + "lensed arcs \u2014 off-arc pixels carry no information about $\\delta\\psi$. The data is therefore masked to an\n", + "arc-tracing region: `al.pc.util.arc_mask_from` thresholds the signal-to-noise map, drops small islands, dilates\n", + "the result and cleans it so every unmasked pixel supports the finite-difference operators built below." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "mask_array = al.pc.util.arc_mask_from(\n", + " np.asarray(dataset.signal_to_noise_map.native), threshold=3.0, ignore_size=10, ext_size=3\n", + ")\n", + "mask = al.Mask2D(mask=mask_array, pixel_scales=dataset.pixel_scales)\n", + "masked_imaging = dataset.apply_mask(mask=mask)\n", + "\n", + "n_data = int(np.count_nonzero(~mask_array))\n", + "print(f\"unmasked data pixels: n_data = {n_data}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Smooth Starting Model__\n", + "\n", + "Potential corrections perturb a smooth starting model \u2014 in a real analysis, the maximum-likelihood result of a\n", + "standard parametric fit. Here we use the true smooth lens (without the subhalo) and the true source.\n", + "\n", + "The source enters the formalism through a `SrcFactory`, which can evaluate the source's brightness and its\n", + "gradients $(\\partial S / \\partial y, \\partial S / \\partial x)$ at arbitrary source-plane positions." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "lens_smooth = al.Galaxy(redshift=0.2, mass=lens_true.mass)\n", + "source_start = al.pc.AnalyticSrcFactory(source_galaxy=source_true)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Dpsi Mesh__\n", + "\n", + "The corrections are defined on a rectangular mesh a factor coarser than the data grid (here factor 2).\n", + "`PairRegularDpsiMesh` builds everything the formalism needs from the data mask alone:\n", + "\n", + " - `mask_dpsi`: the coarse mesh's mask (a coarse pixel is unmasked only if all its data pixels are);\n", + " - `itp_mat`: a sparse [n_data, n_dpsi] bilinear interpolation matrix from mesh to data grid;\n", + " - `Hx_dpsi`, `Hy_dpsi`: sparse first-derivative operators on the mesh's unmasked pixels;\n", + " - `hamiltonian_dpsi`: the mesh Laplacian, used later to convert $\\delta\\psi$ to $\\delta\\kappa$.\n", + "\n", + "Every row of `itp_mat` holds the four bilinear weights of the coarse-mesh box enclosing that data pixel, so its\n", + "rows sum to exactly 1." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "pair = al.pc.PairRegularDpsiMesh(\n", + " mask_array, pixel_scale=dataset.pixel_scales[0], dpsi_factor=2\n", + ")\n", + "\n", + "n_dpsi = int(np.count_nonzero(~pair.mask_dpsi))\n", + "print(f\"unmasked dpsi mesh pixels: n_dpsi = {n_dpsi}\")\n", + "print(f\"itp_mat shape = {pair.itp_mat.shape}, row sums all 1: {np.allclose(pair.itp_mat.sum(axis=1), 1.0)}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Dpsi Gradient Operator__\n", + "\n", + "The correction $\\delta\\psi$ deflects rays by its gradient: $\\delta\\alpha = \\nabla \\delta\\psi$. The sparse operator\n", + "\n", + "$D_\\psi = \\rm{interleave}(\\rm{itp} \\cdot H_x, \\; \\rm{itp} \\cdot H_y)$\n", + "\n", + "of shape [2 n_data, n_dpsi] takes the mesh corrections to their (x, y) gradients at every data pixel, with the\n", + "per-pixel rows interleaved as $(x_0, y_0, x_1, y_1, ...)$." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "dpsi_gradient_matrix = al.pc.util.dpsi_gradient_matrix_from(\n", + " pair.itp_mat, pair.Hx_dpsi, pair.Hy_dpsi\n", + ")\n", + "print(f\"D_psi shape = {dpsi_gradient_matrix.shape}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Source Gradients__\n", + "\n", + "A small extra deflection $\\delta\\alpha$ at an image pixel re-samples the source at a position shifted\n", + "by $-\\delta\\alpha$, changing the observed brightness by $-\\nabla S \\cdot \\delta\\alpha$ to first order. We\n", + "therefore need the source's gradients at the ray-traced (source-plane) position of every image pixel:\n", + "\n", + " - ray-trace the masked grid through the smooth lens;\n", + " - evaluate the source's $(\\partial S/\\partial y, \\partial S/\\partial x)$ there by central differences\n", + " (`SrcFactory.eval_grad`);\n", + " - pack them into the sparse [n_data, 2 n_data] matrix $D_s$, whose row $i$ holds $(\\partial_x S_i, \\partial_y S_i)$\n", + " in the columns matching $D_\\psi$'s interleaved rows." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "traced_grid = masked_imaging.grid.slim - lens_smooth.deflections_yx_2d_from(\n", + " masked_imaging.grid.slim\n", + ")\n", + "source_gradients = source_start.eval_grad(traced_grid[:, 1], traced_grid[:, 0])\n", + "source_gradient_matrix = al.pc.util.source_gradient_matrix_from(source_gradients)\n", + "print(f\"D_s shape = {source_gradient_matrix.shape}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__PSF Blur Matrix__\n", + "\n", + "The correction's brightness response is blurred by the telescope PSF like any other emission. The formalism uses\n", + "the explicit blur matrix $B$ of shape [n_data, n_data]: column $i$ is the PSF kernel centred on pixel $i$,\n", + "restricted to the mask." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "psf_matrix = al.pc.util.psf_matrix_from(\n", + " np.asarray(masked_imaging.psf.kernel.native), np.asarray(masked_imaging.mask)\n", + ")\n", + "print(f\"B shape = {psf_matrix.shape}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Dpsi Mapping Matrix__\n", + "\n", + "Combining the three operators gives the linear response of the observed image to the mesh corrections\n", + "(eq. 8-9 of the potential-correction formalism; Cao et al. 2025):\n", + "\n", + "$\\delta d = - B \\, D_s \\, D_\\psi \\, \\delta\\psi$\n", + "\n", + "The minus sign: a positive potential bump deflects rays outward, sampling the source closer to its centre." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "dpsi_mapping_matrix = np.asarray(-1.0 * psf_matrix @ source_gradient_matrix @ dpsi_gradient_matrix)\n", + "print(f\"dpsi mapping matrix shape = {dpsi_mapping_matrix.shape}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Source Inversion Blocks__\n", + "\n", + "The joint inversion reconstructs the source simultaneously. Its source blocks come from the standard pixelized\n", + "source inversion at the smooth model (see `features/pixelization/likelihood_function.py` for the full\n", + "walkthrough): the PSF-operated mapping matrix $F_{\\rm src}$ and the source regularization matrix $R_s$." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "src_pixelization = al.Pixelization(\n", + " mesh=al.mesh.RectangularUniform(shape=(20, 20)),\n", + " regularization=al.reg.Constant(coefficient=1.0),\n", + ")\n", + "\n", + "source_galaxy = al.Galaxy(redshift=0.6, pixelization=src_pixelization)\n", + "tracer = al.Tracer(galaxies=[lens_smooth, source_galaxy])\n", + "src_fit = al.FitImaging(\n", + " dataset=masked_imaging.apply_over_sampling(\n", + " over_sample_size_lp=4, over_sample_size_pixelization=4\n", + " ),\n", + " tracer=tracer,\n", + " settings=al.Settings(use_positive_only_solver=True, use_border_relocator=True),\n", + ")\n", + "src_mapping_matrix = np.asarray(src_fit.inversion.operated_mapping_matrix)\n", + "src_regularization_matrix = np.asarray(src_fit.inversion.regularization_matrix)\n", + "n_src = src_regularization_matrix.shape[0]\n", + "print(f\"F_src shape = {src_mapping_matrix.shape}, n_src = {n_src}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Joint System__\n", + "\n", + "The joint linear system stacks the two blocks:\n", + "\n", + " - mapping matrix $M = [\\, F_{\\rm src} \\; | \\; -B D_s D_\\psi \\,]$ of shape [n_data, n_src + n_dpsi];\n", + " - block-diagonal regularization $R = \\rm{diag}(R_s, R_\\psi)$, where $R_\\psi$ is the `MaternKernel` scheme built\n", + " from the dpsi mesh's unmasked pixel positions through the `DpsiLinearObj` adapter (mask-based schemes like\n", + " `CurvatureMask` plug into the same adapter, which exposes both the mesh's `mask` and its pixel positions).\n", + "\n", + "Because the two blocks share one linear solve, the source-vs-corrections covariance is fully accounted for \u2014 the\n", + "corrections cannot silently absorb source structure the data attributes to the source, and vice versa." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "dpsi_regularization = al.reg.MaternKernel(coefficient=2000.0, scale=4.0, nu=2.5)\n", + "dpsi_points = np.vstack([pair.ygrid_dpsi_1d, pair.xgrid_dpsi_1d]).T\n", + "dpsi_linear_obj = al.pc.DpsiLinearObj(mask=pair.mask_dpsi, points=dpsi_points)\n", + "dpsi_regularization_matrix = dpsi_regularization.regularization_matrix_from(\n", + " linear_obj=dpsi_linear_obj\n", + ")\n", + "\n", + "mapping_matrix = np.hstack([src_mapping_matrix, dpsi_mapping_matrix])\n", + "regularization_matrix = np.asarray(\n", + " block_diag([src_regularization_matrix, dpsi_regularization_matrix]).toarray()\n", + ")\n", + "print(f\"joint mapping matrix shape = {mapping_matrix.shape}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Solve__\n", + "\n", + "With diagonal noise covariance $C^{-1} = \\rm{diag}(1/\\sigma_i^2)$, the maximum-evidence solution of the joint\n", + "system is the standard regularized normal equations (Warren & Dye 2003; Suyu et al. 2006):\n", + "\n", + "$(M^T C^{-1} M + R) \\, x = M^T C^{-1} d$\n", + "\n", + "whose solution vector splits into the reconstructed source (first n_src entries) and the corrections\n", + "(last n_dpsi entries)." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "data = np.asarray(masked_imaging.data)\n", + "noise = np.asarray(masked_imaging.noise_map)\n", + "inv_variance = 1.0 / noise**2\n", + "\n", + "curvature_matrix = mapping_matrix.T @ (mapping_matrix * inv_variance[:, None])\n", + "data_vector = mapping_matrix.T @ (inv_variance * data)\n", + "curvature_reg_matrix = curvature_matrix + regularization_matrix\n", + "\n", + "solution = np.linalg.solve(curvature_reg_matrix, data_vector)\n", + "source_solution = solution[:n_src]\n", + "dpsi_solution = solution[n_src:]\n", + "model_image = mapping_matrix @ solution" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Dkappa Map__\n", + "\n", + "The corrections' physical meaning is clearest as a convergence correction,\n", + "$\\delta\\kappa = \\frac{1}{2}\\nabla^2 \\delta\\psi$, via the mesh Laplacian. A dark subhalo missing from the smooth\n", + "model appears as a positive $\\delta\\kappa$ peak at its position." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "dkappa = pair.hamiltonian_dpsi @ dpsi_solution\n", + "peak = int(np.argmax(dkappa))\n", + "print(\n", + " f\"dkappa peak at (y, x) = ({pair.ygrid_dpsi_1d[peak]:.2f}, {pair.xgrid_dpsi_1d[peak]:.2f})\"\n", + " f\" \u2014 true subhalo at (1.41, 0.00)\"\n", + ")\n", + "\n", + "fig, axes = plt.subplots(1, 2, figsize=(10, 4))\n", + "al.pc.visualize.imshow_masked_data(\n", + " dpsi_solution, pair.mask_dpsi, ax=axes[0], origin=\"upper\", extent=pair.data_bound\n", + ")\n", + "axes[0].set_title(\"dpsi\")\n", + "al.pc.visualize.imshow_masked_data(\n", + " dkappa, pair.mask_dpsi, ax=axes[1], origin=\"upper\", extent=pair.data_bound\n", + ")\n", + "axes[1].set_title(\"dkappa\")\n", + "plt.tight_layout()\n", + "plt.show()\n", + "plt.close()" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Evidence Terms__\n", + "\n", + "The Bayesian evidence of the joint inversion has five terms (Suyu et al. 2006 eq. 19; Cao et al. 2025):\n", + "\n", + " 1. the noise normalization $-\\frac{1}{2}\\sum_i \\log(2\\pi\\sigma_i^2)$;\n", + " 2. the Occam term $-\\frac{1}{2}\\log\\det(M^T C^{-1} M + R)$, penalising flexible models;\n", + " 3. the regularization normalizations $+\\frac{1}{2}[\\log\\det R_s + \\log\\det R_\\psi]$;\n", + " 4. the regularization penalty of the solution $-\\frac{1}{2} x^T R x$;\n", + " 5. the $\\chi^2$ of the fit $-\\frac{1}{2}\\sum_i (d_i - m_i)^2/\\sigma_i^2$.\n", + "\n", + "Terms 2-4 are what allow the evidence to set both regularization strengths objectively: stronger regularization\n", + "lowers the Occam term but raises the penalty and chi-squared, and the evidence peaks at the balance." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "noise_term = -0.5 * float(np.sum(np.log(2.0 * np.pi * noise**2)))\n", + "occam_term = -0.5 * float(np.linalg.slogdet(curvature_reg_matrix)[1])\n", + "reg_norm_term = 0.5 * (\n", + " float(np.linalg.slogdet(src_regularization_matrix)[1])\n", + " + float(np.linalg.slogdet(dpsi_regularization_matrix)[1])\n", + ")\n", + "reg_penalty_term = -0.5 * float(solution @ regularization_matrix @ solution)\n", + "chi_squared_term = -0.5 * float(np.sum(((data - model_image) / noise) ** 2))\n", + "\n", + "log_evidence = (\n", + " noise_term + occam_term + reg_norm_term + reg_penalty_term + chi_squared_term\n", + ")\n", + "\n", + "print(f\"noise term = {noise_term:.6e}\")\n", + "print(f\"occam term = {occam_term:.6e}\")\n", + "print(f\"reg normalizations = {reg_norm_term:.6e}\")\n", + "print(f\"reg penalty = {reg_penalty_term:.6e}\")\n", + "print(f\"chi squared term = {chi_squared_term:.6e}\")\n", + "print(f\"log evidence = {log_evidence:.8e}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Verification__\n", + "\n", + "`al.pc.FitDpsiSrcImaging` performs exactly the steps above. Its `log_evidence` must equal our explicit\n", + "calculation (we preload the source blocks so both use the identical source inversion)." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "fit = al.pc.FitDpsiSrcImaging(\n", + " masked_imaging=masked_imaging,\n", + " lens_start=lens_smooth,\n", + " source_start=source_start,\n", + " dpsi_pixelization=al.pc.DpsiPixelization(\n", + " mesh=al.pc.RegularDpsiMesh(factor=2), regularization=dpsi_regularization\n", + " ),\n", + " src_pixelization=src_pixelization,\n", + " preloads={\n", + " \"src_map_mat\": src_mapping_matrix,\n", + " \"src_reg_mat\": src_regularization_matrix,\n", + " },\n", + ")\n", + "\n", + "print(f\"al.pc.FitDpsiSrcImaging.log_evidence = {fit.log_evidence:.8e}\")\n", + "assert np.isclose(fit.log_evidence, log_evidence, rtol=1e-8)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__JAX / xp__\n", + "\n", + "Every dense step above is also available through `al.pc.dense_util`, whose kernels follow the PyAuto `xp`\n", + "convention: written once with `xp=np` (what we just did by hand), and runnable under `xp=jax.numpy` for\n", + "jit-compiled, accelerator-ready execution \u2014 see\n", + "`autolens_workspace_test/scripts/jax_likelihood_functions/imaging/potential_correction.py` for the numerical\n", + "agreement checks.\n", + "\n", + " from autolens.potential_correction import dense_util\n", + "\n", + " result = dense_util.log_evidence_joint_dense_from(\n", + " data, noise, mapping_matrix, src_regularization_matrix, dpsi_regularization_matrix, xp=np, # or jax.numpy\n", + " )\n", + "\n", + "__Wrap Up__\n", + "\n", + "We built the potential-correction evidence from its raw ingredients: the dpsi mesh and its interpolation to the\n", + "data grid, the gradient operators, the PSF blur matrix, the joint block system and its five evidence terms \u2014\n", + "and verified `al.pc.FitDpsiSrcImaging` reproduces the number exactly.\n", + "\n", + "In a science analysis the regularization hyper-parameters are sampled with a non-linear search through\n", + "`al.pc.DpsiSrcInvAnalysis` (one-shot) or refined with the iterative Levenberg-Marquardt engine\n", + "`al.pc.IterFitDpsiSrcImaging` \u2014 see `guides/advanced/potential_correction.py`.\n", + "\n", + "If you use this functionality, please cite Cao et al. 2025\n", + "(https://github.com/caoxiaoyue/potential_correction_paper) alongside PyAutoLens." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [], + "outputs": [], + "execution_count": null + } + ], + "metadata": { + "anaconda-cloud": {}, + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.1" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} \ No newline at end of file diff --git a/notebooks/interferometer/features/README.md b/notebooks/interferometer/features/README.md index f518afcbd..bccb3d6a1 100644 --- a/notebooks/interferometer/features/README.md +++ b/notebooks/interferometer/features/README.md @@ -12,6 +12,7 @@ The following example scripts illustrating lens modeling where: - `pixelization`: The source is reconstructed using an adaptive rectangular or Delaunay mesh - `extra_galaxies`: Modeling which account for the light and mass of extra nearby galaxies. - `subhalo`: Fitting lens models for dark matter subhalo detection and sensitivity mapping. +- `potential_correction`: Gravitational imaging — pixelized corrections to the lensing potential reconstructed jointly with the source, via the visibility-count-independent sparse-operator route. # Notes diff --git a/notebooks/interferometer/features/potential_correction/likelihood_function.ipynb b/notebooks/interferometer/features/potential_correction/likelihood_function.ipynb new file mode 100644 index 000000000..f159c6ed9 --- /dev/null +++ b/notebooks/interferometer/features/potential_correction/likelihood_function.ipynb @@ -0,0 +1,676 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Log Likelihood Function: Potential Correction (Interferometer)__\n", + "\n", + "This script provides a step-by-step guide of the **PyAutoLens** visibility-space potential-correction\n", + "`log_likelihood_function` (the Bayesian evidence of `al.pc.FitDpsiSrcInterferometer`), which jointly inverts\n", + "`Interferometer` data for a pixelized source and pixelized corrections $\\delta\\psi$ to the lensing potential.\n", + "\n", + "This script has the following aims:\n", + "\n", + " - To provide a resource that authors can include in papers using **PyAutoLens**, so that readers can understand the\n", + " likelihood function (including references to the previous literature from which it is defined) without having to\n", + " write large quantities of text and equations.\n", + "\n", + " - To make visibility-space gravitational imaging less of a \"black-box\" to users: every operator of the correction\n", + " formalism ($D_s$, $D_\\psi$, the joint response and both curvature routes) is built explicitly as a numpy array below.\n", + "\n", + "If you use the potential-correction functionality in your research, please cite Cao et al. 2025, from whose\n", + "`potential_correction` package (https://github.com/caoxiaoyue/lensing_potential_correction) the implementation is\n", + "ported; citation materials are provided at https://github.com/caoxiaoyue/potential_correction_paper. The method\n", + "builds on Koopmans 2005, Suyu et al. 2009 and Vegetti & Koopmans 2009; the visibility-space formulation is the\n", + "methodology behind the JVAS B1938+666 detections of Powell et al. 2025 (Nature Astronomy 9, 1714) and\n", + "Vegetti et al. 2026.\n", + "\n", + "__Contents__\n", + "\n", + "- **Simplifications:** The choices made to keep this walkthrough small and explicit.\n", + "- **Prerequisites:** The likelihood functions this one builds on.\n", + "- **Dataset:** Simulate the interferometer dataset fitted: a lens whose true mass contains a dark subhalo.\n", + "- **Source Inversion Blocks:** The standard pixelized source inversion supplying the source blocks.\n", + "- **Dpsi Mesh:** The arc-restricted coarse mesh the corrections are defined on.\n", + "- **Correction Response:** The real-space response G = -D_s D_psi of the image to the corrections.\n", + "- **Joint Response:** The joint real-space response A = [f | G] \u2014 no PSF matrix: the measurement operator here\n", + " is the non-uniform Fourier transform.\n", + "- **Dense Route:** The visibility-space normal equations built from the explicitly transformed response.\n", + "- **Sparse (w-tilde) Route:** The same curvature and data vector without ever forming the transformed response.\n", + "- **Solve:** Solving the joint linear system for the source and the corrections.\n", + "- **Dkappa Map:** Converting the solved corrections into a convergence-correction map.\n", + "- **Evidence Terms:** The five terms of the Bayesian evidence, computed explicitly.\n", + "- **Verification:** The same number from `al.pc.FitDpsiSrcInterferometer.log_evidence`, both routes.\n", + "- **JAX / xp:** The xp-API kernels of the same computation.\n", + "- **Wrap Up:** Summary and next steps.\n", + "\n", + "__Simplifications__\n", + "\n", + "To keep every array small enough to inspect, this example uses a low-resolution simulated dataset with only 1000\n", + "visibilities \u2014 small enough that the dense reference route (which materializes the [2 n_vis, n_par] transformed\n", + "response) fits in memory alongside the sparse route. The source is reconstructed on a `RectangularUniform` mesh\n", + "(uniform source pixels are simpler to reason about than adaptive meshes) and the corrections are regularized with\n", + "the `MaternKernel` scheme (`nu=2.5`) used in science analyses. The smooth starting model and source are the truths\n", + "used in the simulation, isolating the correction formalism itself.\n", + "\n", + "__Prerequisites__\n", + "\n", + "- `interferometer/features/pixelization/likelihood_function.ipynb` \u2014 the visibility-space pixelized-source\n", + " evidence, whose terms reappear here.\n", + "- `imaging/features/potential_correction/likelihood_function.ipynb` \u2014 the imaging analogue of this walkthrough,\n", + " where the measurement operator is the PSF-blur matrix instead of the Fourier transform.\n", + "- `interferometer/features/potential_correction/start_here.ipynb` \u2014 the user-facing example of this feature." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Google Colab Setup__\n", + "\n", + "This cell sets up the environment when the notebook is run on Google Colab: it installs the\n", + "required PyAuto packages, clones the workspace (configuration files and example datasets) and\n", + "points the configuration at it. If you are running the notebook elsewhere (e.g. locally via\n", + "your own installation) it does nothing, and you can run it safely.\n", + "\n", + "Colab tip: model-fits run much faster on a GPU \u2014 enable one via \"Runtime\" -> \"Change runtime\n", + "type\" -> \"Hardware accelerator\" before running the notebook." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "try:\n", + " import google.colab\n", + " import subprocess\n", + " import sys\n", + "\n", + " subprocess.check_call(\n", + " [sys.executable, \"-m\", \"pip\", \"install\", \"autoconf\", \"--no-deps\"]\n", + " )\n", + "except ImportError:\n", + " pass\n", + "\n", + "from autoconf import setup_colab\n", + "\n", + "setup_colab.setup(\"autolens\")" + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "\n", + "from autoconf import jax_wrapper # Sets JAX environment before other imports\n", + "\n", + "from autoconf import setup_notebook; setup_notebook()\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from scipy.sparse import block_diag\n", + "\n", + "import autoarray as aa\n", + "import autolens as al\n", + "import autolens.plot as aplt" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Dataset__\n", + "\n", + "We simulate the dataset in-memory (seeded, so this script is fully reproducible): random uv coverage, an\n", + "`Isothermal` lens whose true mass also contains a $10^{10} M_\\odot$ NFW subhalo on the Einstein ring, and a\n", + "compact double-Gaussian source. The subhalo is what the corrections will recover.\n", + "\n", + "Unlike imaging, the real-space mask must stay a **filled circle**: it defines the rectangular extent of the FFTs\n", + "inside the sparse operator." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "rng = np.random.default_rng(1)\n", + "uv_wavelengths = rng.uniform(-3.0e5, 3.0e5, size=(1000, 2))\n", + "\n", + "real_space_mask = al.Mask2D.circular(shape_native=(64, 64), pixel_scales=0.1, radius=2.6)\n", + "grid = al.Grid2D.from_mask(mask=real_space_mask)\n", + "\n", + "simulator = al.SimulatorInterferometer(\n", + " uv_wavelengths=uv_wavelengths,\n", + " exposure_time=300.0,\n", + " noise_sigma=0.02,\n", + " noise_seed=1,\n", + ")\n", + "\n", + "lens_true = al.Galaxy(\n", + " redshift=0.2,\n", + " mass=al.mp.Isothermal(\n", + " centre=(0.0, 0.0),\n", + " einstein_radius=1.4,\n", + " ell_comps=al.convert.ell_comps_from(axis_ratio=0.9, angle=0.0),\n", + " ),\n", + " subhalo=al.mp.NFWMCRLudlowSph(\n", + " centre=(1.41, 0.0), mass_at_200=1.0e10, redshift_object=0.2, redshift_source=0.6\n", + " ),\n", + ")\n", + "source_true = al.Galaxy(\n", + " redshift=0.6,\n", + " bulge0=al.lp.Gaussian(\n", + " centre=(0.0, 0.0),\n", + " ell_comps=al.convert.ell_comps_from(axis_ratio=0.6, angle=45.0),\n", + " intensity=5.0,\n", + " sigma=0.15,\n", + " ),\n", + " bulge1=al.lp.Gaussian(\n", + " centre=(0.0, 0.4),\n", + " ell_comps=al.convert.ell_comps_from(axis_ratio=0.4, angle=135.0),\n", + " intensity=3.0,\n", + " sigma=0.1,\n", + " ),\n", + ")\n", + "\n", + "dataset = simulator.via_tracer_from(\n", + " tracer=al.Tracer(galaxies=[lens_true, source_true]), grid=grid\n", + ")\n", + "dataset = al.Interferometer(\n", + " data=dataset.data,\n", + " noise_map=dataset.noise_map,\n", + " uv_wavelengths=uv_wavelengths,\n", + " real_space_mask=real_space_mask,\n", + ")\n", + "dataset = dataset.apply_sparse_operator()\n", + "\n", + "n_full = int(np.count_nonzero(~np.asarray(real_space_mask)))\n", + "print(f\"visibilities: {uv_wavelengths.shape[0]}, real-space mask pixels: n_full = {n_full}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Smooth Starting Model__\n", + "\n", + "Potential corrections perturb a smooth starting model \u2014 in a real analysis, the maximum-likelihood result of a\n", + "standard parametric fit. Here we use the true smooth lens (without the subhalo) and the true source.\n", + "\n", + "The source enters the formalism through a `SrcFactory`, which can evaluate the source's brightness and its\n", + "gradients $(\\partial S / \\partial y, \\partial S / \\partial x)$ at arbitrary source-plane positions." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "lens_smooth = al.Galaxy(redshift=0.2, mass=lens_true.mass)\n", + "source_start = al.pc.AnalyticSrcFactory(source_galaxy=source_true)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Source Inversion Blocks__\n", + "\n", + "The joint inversion reconstructs the source simultaneously. Its source blocks come from the standard pixelized\n", + "visibility-space source inversion at the smooth model (see\n", + "`interferometer/features/pixelization/likelihood_function.py` for the full walkthrough): the mapper's real-space\n", + "mapping matrix $f$ of shape [n_full, n_src] and the source regularization matrix $R_s$.\n", + "\n", + "Note that $f$ is a **real-space** matrix \u2014 unlike imaging, where the source block is PSF-convolved before entering\n", + "the joint system, here the measurement operator (the Fourier transform) is applied to the whole joint response at\n", + "once, further below." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "src_pixelization = al.Pixelization(\n", + " mesh=al.mesh.RectangularUniform(shape=(20, 20)),\n", + " regularization=al.reg.Constant(coefficient=3.8),\n", + ")\n", + "source_galaxy = al.Galaxy(redshift=1.0, pixelization=src_pixelization)\n", + "src_fit = al.FitInterferometer(\n", + " dataset=dataset,\n", + " tracer=al.Tracer(galaxies=[lens_smooth, source_galaxy]),\n", + " settings=aa.Settings(use_positive_only_solver=True, use_border_relocator=True),\n", + ")\n", + "mapper = src_fit.inversion.linear_obj_list[0]\n", + "src_mapping_matrix = np.asarray(mapper.mapping_matrix)\n", + "src_regularization_matrix = np.asarray(src_fit.inversion.regularization_matrix)\n", + "n_src = src_regularization_matrix.shape[0]\n", + "print(f\"f shape = {src_mapping_matrix.shape}, n_src = {n_src}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Dpsi Mesh__\n", + "\n", + "The corrections are only constrained where the lensed arcs are, so their mesh is restricted to an arc-tracing\n", + "sub-mask of the real-space mask (built here from the smooth-model image geometry). `PairRegularDpsiMesh` then\n", + "builds everything the formalism needs on a mesh a factor coarser than the data grid:\n", + "\n", + " - `mask_data` / `mask_dpsi`: the cleaned arc mask on the data grid and its coarse-mesh counterpart;\n", + " - `itp_mat`: a sparse bilinear interpolation matrix from mesh to data grid (rows sum to exactly 1);\n", + " - `Hx_dpsi`, `Hy_dpsi`: sparse first-derivative operators on the mesh's unmasked pixels;\n", + " - `hamiltonian_dpsi`: the mesh Laplacian, used later to convert $\\delta\\psi$ to $\\delta\\kappa$." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "tracer_smooth = al.Tracer(galaxies=[lens_smooth, source_true])\n", + "arc_image = np.asarray(tracer_smooth.image_2d_from(grid=grid).native)\n", + "arc_mask = al.pc.util.arc_mask_from(\n", + " arc_image / (0.05 * arc_image.max()), threshold=3.0, ignore_size=10, ext_size=3\n", + ")\n", + "dpsi_mask = ~((~arc_mask) & (~np.asarray(real_space_mask)))\n", + "\n", + "pair = al.pc.PairRegularDpsiMesh(dpsi_mask, pixel_scale=0.1, dpsi_factor=2)\n", + "n_sub = int(np.count_nonzero(~pair.mask_data))\n", + "n_dpsi = int(np.count_nonzero(~pair.mask_dpsi))\n", + "print(f\"arc pixels: n_sub = {n_sub}, dpsi mesh pixels: n_dpsi = {n_dpsi}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Correction Response__\n", + "\n", + "As in the imaging walkthrough, a correction $\\delta\\psi$ perturbs the observed image via the source's brightness\n", + "gradients at the ray-traced position of every (arc) image pixel:\n", + "\n", + "$\\delta I = - D_s \\, D_\\psi \\, \\delta\\psi \\equiv G \\, \\delta\\psi$\n", + "\n", + " - $D_\\psi = \\rm{interleave}(\\rm{itp} \\cdot H_x, \\; \\rm{itp} \\cdot H_y)$, of shape [2 n_sub, n_dpsi], takes the\n", + " mesh corrections to their (x, y) gradients at every arc pixel;\n", + " - $D_s$, of shape [n_sub, 2 n_sub], holds the source's gradients at the traced positions;\n", + " - the minus sign: a positive potential bump deflects rays outward, sampling the source closer to its centre.\n", + "\n", + "The arc rows are then scattered into the full real-space row space (zero response off the arcs), giving the\n", + "real-space correction response $G$ of shape [n_full, n_dpsi]. There is **no PSF-blur matrix**: blurring by the\n", + "instrument response happens in visibility space, through the Fourier transform below." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "dpsi_gradient_matrix = al.pc.util.dpsi_gradient_matrix_from(\n", + " pair.itp_mat, pair.Hx_dpsi, pair.Hy_dpsi\n", + ")\n", + "\n", + "traced_grid = np.asarray(\n", + " dataset.grid.slim - lens_smooth.deflections_yx_2d_from(dataset.grid.slim)\n", + ")\n", + "full_index = np.full(np.asarray(real_space_mask).shape, -1, dtype=int)\n", + "full_index[~np.asarray(real_space_mask)] = np.arange(n_full)\n", + "rows_in_full = full_index[~pair.mask_data]\n", + "\n", + "source_gradients = source_start.eval_grad(\n", + " traced_grid[rows_in_full, 1], traced_grid[rows_in_full, 0]\n", + ")\n", + "source_gradient_matrix = al.pc.util.source_gradient_matrix_from(source_gradients)\n", + "\n", + "G_sub = np.asarray((-1.0 * source_gradient_matrix @ dpsi_gradient_matrix).todense())\n", + "G = np.zeros((n_full, n_dpsi))\n", + "G[rows_in_full] = G_sub\n", + "print(f\"G shape = {G.shape}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Joint Response__\n", + "\n", + "The joint real-space response stacks the two blocks, $A = [\\, f \\; | \\; G \\,]$, alongside the block-diagonal\n", + "regularization $R = \\rm{diag}(R_s, R_\\psi)$, where $R_\\psi$ is the `MaternKernel` scheme built from the dpsi\n", + "mesh's unmasked pixel positions through the `DpsiLinearObj` adapter.\n", + "\n", + "Because the two blocks share one linear solve, the source-vs-corrections covariance is fully accounted for \u2014 the\n", + "corrections cannot silently absorb structure the data attributes to the source, and vice versa." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "dpsi_regularization = al.reg.MaternKernel(coefficient=2000.0, scale=4.0, nu=2.5)\n", + "dpsi_points = np.vstack([pair.ygrid_dpsi_1d, pair.xgrid_dpsi_1d]).T\n", + "dpsi_linear_obj = al.pc.DpsiLinearObj(mask=pair.mask_dpsi, points=dpsi_points)\n", + "dpsi_regularization_matrix = np.asarray(\n", + " dpsi_regularization.regularization_matrix_from(linear_obj=dpsi_linear_obj)\n", + ")\n", + "\n", + "A = np.hstack([src_mapping_matrix, G])\n", + "regularization_matrix = np.asarray(\n", + " block_diag([src_regularization_matrix, dpsi_regularization_matrix]).toarray()\n", + ")\n", + "print(f\"joint real-space response A shape = {A.shape}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Dense Route__\n", + "\n", + "The measurement operator is the non-uniform Fourier transform $T$, taking any real-space image to model\n", + "visibilities. The dense reference route materializes the transformed joint response $T(A)$, stacks its real and\n", + "imaginary parts row-wise into $M$ of shape [2 n_vis, n_par], and forms the standard regularized normal equations\n", + "(Warren & Dye 2003; Suyu et al. 2006) with diagonal noise covariance $C^{-1}$:\n", + "\n", + "curvature $F = M^T C^{-1} M$, data vector $D = M^T C^{-1} d$.\n", + "\n", + "This is exact but scales with the visibility count \u2014 fine at 1000 visibilities, impossible at the $10^5$-$10^8$\n", + "of real ALMA/VLBI datasets." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "transformed = dataset.transformer.transform_mapping_matrix(A)\n", + "M = np.vstack([np.real(transformed), np.imag(transformed)])\n", + "\n", + "data = np.asarray(dataset.data)\n", + "noise = np.asarray(dataset.noise_map)\n", + "stacked_data = np.concatenate([data.real, data.imag])\n", + "stacked_inv_variance = 1.0 / np.concatenate([noise.real, noise.imag]) ** 2\n", + "\n", + "curvature_matrix = M.T @ (M * stacked_inv_variance[:, None])\n", + "data_vector = M.T @ (stacked_inv_variance * stacked_data)\n", + "print(f\"dense route: M shape = {M.shape}, curvature shape = {curvature_matrix.shape}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Sparse (w-tilde) Route__\n", + "\n", + "The production route computes the identical $F$ and $D$ without ever forming $T(A)$ (Powell et al. 2021's\n", + "visibility-space w-tilde formalism):\n", + "\n", + " - the curvature is $F = A^T \\, (T^H C^{-1} T) \\, A$, where the operator $T^H C^{-1} T$ is a **convolution** in\n", + " real space \u2014 applied by FFTs on the mask's rectangular extent, so the cost depends on the number of mask\n", + " pixels, not the number of visibilities;\n", + " - the data vector is $D = A^T \\tilde{d}$, where $\\tilde{d} = T^H C^{-1} d$ is the **dirty image** of the\n", + " visibilities \u2014 computed once when `apply_sparse_operator` is called.\n", + "\n", + "Both are available from the dataset's sparse operator; we verify they equal the dense-route matrices." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "data_vector_sparse = A.T @ np.asarray(dataset.sparse_operator.dirty_image)\n", + "print(\n", + " f\"sparse data vector == dense data vector: \"\n", + " f\"{np.allclose(data_vector_sparse, data_vector, rtol=1e-6)}\"\n", + ")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The sparse curvature is assembled from the COO triplets of $A$ (the machinery `al.pc.FitDpsiSrcInterferometer`\n", + "uses internally); rather than repeat that plumbing here, we take the fit's sparse-route curvature below and verify\n", + "it against our dense $F$ at the end. For this walkthrough we continue with the dense-route matrices, which we have\n", + "just shown are the same thing.\n", + "\n", + "__Solve__\n", + "\n", + "The maximum-evidence solution of the joint system:\n", + "\n", + "$(F + R) \\, x = D$\n", + "\n", + "splits into the reconstructed source (first n_src entries) and the corrections (last n_dpsi entries)." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "curvature_reg_matrix = curvature_matrix + regularization_matrix\n", + "solution = np.linalg.solve(curvature_reg_matrix, data_vector)\n", + "source_solution = solution[:n_src]\n", + "dpsi_solution = solution[n_src:]" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Dkappa Map__\n", + "\n", + "The corrections' physical meaning is clearest as a convergence correction,\n", + "$\\delta\\kappa = \\frac{1}{2}\\nabla^2 \\delta\\psi$, via the mesh Laplacian. A dark subhalo missing from the smooth\n", + "model appears as a positive $\\delta\\kappa$ peak at its position." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "dkappa = pair.hamiltonian_dpsi @ dpsi_solution\n", + "peak = int(np.argmax(dkappa))\n", + "print(\n", + " f\"dkappa peak at (y, x) = ({pair.ygrid_dpsi_1d[peak]:.2f}, {pair.xgrid_dpsi_1d[peak]:.2f})\"\n", + " f\" \u2014 true subhalo at (1.41, 0.00)\"\n", + ")\n", + "\n", + "fig, axes = plt.subplots(1, 2, figsize=(10, 4))\n", + "al.pc.visualize.imshow_masked_data(\n", + " dpsi_solution, pair.mask_dpsi, ax=axes[0], origin=\"upper\", extent=pair.data_bound\n", + ")\n", + "axes[0].set_title(\"dpsi\")\n", + "al.pc.visualize.imshow_masked_data(\n", + " dkappa, pair.mask_dpsi, ax=axes[1], origin=\"upper\", extent=pair.data_bound\n", + ")\n", + "axes[1].set_title(\"dkappa\")\n", + "plt.tight_layout()\n", + "plt.show()\n", + "plt.close()" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Evidence Terms__\n", + "\n", + "The Bayesian evidence of the joint inversion has five terms (Suyu et al. 2006 eq. 19; Cao et al. 2025), with the\n", + "visibility-space specifics that the noise normalization and $\\chi^2$ run over the real **and** imaginary parts:\n", + "\n", + " 1. the noise normalization $-\\frac{1}{2}\\sum_i [\\log(2\\pi\\sigma_{R,i}^2) + \\log(2\\pi\\sigma_{I,i}^2)]$;\n", + " 2. the Occam term $-\\frac{1}{2}\\log\\det(F + R)$, penalising flexible models;\n", + " 3. the regularization normalizations $+\\frac{1}{2}[\\log\\det R_s + \\log\\det R_\\psi]$;\n", + " 4. the regularization penalty of the solution $-\\frac{1}{2} x^T R x$;\n", + " 5. the $\\chi^2$ of the fit, $-\\frac{1}{2}\\sum_i |d_i - m_i|^2/\\sigma_i^2$, where the model visibilities\n", + " $m = T(A x)$ cost one forward transform of the reconstructed real-space image \u2014 the only place a transform of\n", + " a model is ever needed in the sparse route.\n", + "\n", + "Terms 2-4 are what allow the evidence to set both regularization strengths objectively: stronger regularization\n", + "lowers the Occam term but raises the penalty and chi-squared, and the evidence peaks at the balance." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "noise_term = -0.5 * float(aa.util.fit.noise_normalization_complex_from(noise_map=noise))\n", + "occam_term = -0.5 * float(np.linalg.slogdet(curvature_reg_matrix)[1])\n", + "reg_norm_term = 0.5 * (\n", + " float(np.linalg.slogdet(src_regularization_matrix)[1])\n", + " + float(np.linalg.slogdet(dpsi_regularization_matrix)[1])\n", + ")\n", + "reg_penalty_term = -0.5 * float(solution @ regularization_matrix @ solution)\n", + "\n", + "model_image = aa.Array2D(values=A @ solution, mask=real_space_mask)\n", + "model_visibilities = np.asarray(dataset.transformer.visibilities_from(image=model_image))\n", + "residual = data - model_visibilities\n", + "chi_squared_term = -0.5 * float(\n", + " np.sum((residual.real / noise.real) ** 2) + np.sum((residual.imag / noise.imag) ** 2)\n", + ")\n", + "\n", + "log_evidence = (\n", + " noise_term + occam_term + reg_norm_term + reg_penalty_term + chi_squared_term\n", + ")\n", + "\n", + "print(f\"noise term = {noise_term:.6e}\")\n", + "print(f\"occam term = {occam_term:.6e}\")\n", + "print(f\"reg normalizations = {reg_norm_term:.6e}\")\n", + "print(f\"reg penalty = {reg_penalty_term:.6e}\")\n", + "print(f\"chi squared term = {chi_squared_term:.6e}\")\n", + "print(f\"log evidence = {log_evidence:.8e}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Verification__\n", + "\n", + "`al.pc.FitDpsiSrcInterferometer` performs exactly the steps above, through either route. Its `log_evidence` must\n", + "equal our explicit calculation for both (we preload the source blocks so all three use the identical source\n", + "inversion), and its sparse-route curvature must equal our dense-route $F$." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "dpsi_pixelization = al.pc.DpsiPixelization(\n", + " mesh=al.pc.RegularDpsiMesh(factor=2), regularization=dpsi_regularization\n", + ")\n", + "preloads = {\"src_mapper\": mapper, \"src_reg_mat\": src_regularization_matrix}\n", + "\n", + "fit_sparse = al.pc.FitDpsiSrcInterferometer(\n", + " dataset=dataset,\n", + " lens_start=lens_smooth,\n", + " source_start=source_start,\n", + " dpsi_pixelization=dpsi_pixelization,\n", + " src_pixelization=src_pixelization,\n", + " dpsi_mask=dpsi_mask,\n", + " use_sparse_operator=True,\n", + " preloads=preloads,\n", + ")\n", + "fit_dense = al.pc.FitDpsiSrcInterferometer(\n", + " dataset=dataset,\n", + " lens_start=lens_smooth,\n", + " source_start=source_start,\n", + " dpsi_pixelization=dpsi_pixelization,\n", + " src_pixelization=src_pixelization,\n", + " dpsi_mask=dpsi_mask,\n", + " use_sparse_operator=False,\n", + " preloads=preloads,\n", + ")\n", + "\n", + "print(f\"sparse-route log_evidence = {fit_sparse.log_evidence:.8e}\")\n", + "print(f\"dense-route log_evidence = {fit_dense.log_evidence:.8e}\")\n", + "\n", + "assert np.allclose(np.asarray(fit_sparse.curvature_matrix), curvature_matrix, rtol=1e-6)\n", + "assert np.isclose(fit_dense.log_evidence, log_evidence, rtol=1e-8)\n", + "assert np.isclose(fit_sparse.log_evidence, log_evidence, rtol=1e-6)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__JAX / xp__\n", + "\n", + "The iterative refinement engine (`al.pc.IterFitDpsiSrcInterferometer`, which re-ray-traces through the corrected\n", + "lens and re-solves) runs its Levenberg-Marquardt kernels through the PyAuto `xp` convention: written once with\n", + "`xp=np` (what we just did by hand), and runnable under `xp=jax.numpy` for jit-compiled, accelerator-ready\n", + "execution. The sparse operator's FFT machinery is JAX-accelerated in all cases \u2014 see\n", + "`autolens_workspace_test/scripts/jax_likelihood_functions/interferometer/potential_correction.py` for the\n", + "numerical agreement checks.\n", + "\n", + "__Wrap Up__\n", + "\n", + "We built the visibility-space potential-correction evidence from its raw ingredients: the arc-restricted dpsi\n", + "mesh, the gradient operators, the joint real-space response, the dense and sparse (w-tilde) routes to the same\n", + "normal equations, and the five evidence terms \u2014 and verified `al.pc.FitDpsiSrcInterferometer` reproduces the\n", + "number exactly through both routes.\n", + "\n", + "In a science analysis the regularization hyper-parameters are sampled with a non-linear search through\n", + "`al.pc.DpsiSrcInvInterferometerAnalysis` (one-shot) or refined with the iterative engine \u2014 see `start_here.py`\n", + "in this folder for the certified end-to-end recipe.\n", + "\n", + "If you use this functionality, please cite Cao et al. 2025\n", + "(https://github.com/caoxiaoyue/potential_correction_paper) alongside PyAutoLens." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [], + "outputs": [], + "execution_count": null + } + ], + "metadata": { + "anaconda-cloud": {}, + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.1" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} \ No newline at end of file diff --git a/notebooks/interferometer/features/potential_correction/start_here.ipynb b/notebooks/interferometer/features/potential_correction/start_here.ipynb new file mode 100644 index 000000000..c51987d12 --- /dev/null +++ b/notebooks/interferometer/features/potential_correction/start_here.ipynb @@ -0,0 +1,530 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Feature: Potential Correction (Gravitational Imaging)\n", + "=====================================================\n", + "\n", + "This example performs gravitational imaging on interferometer data: pixelized corrections $\\delta\\psi$ to the\n", + "lensing potential are reconstructed jointly with the pixelized source, revealing mass structure (e.g. dark\n", + "subhaloes) the smooth lens model omits, as a convergence-correction map $\\delta\\kappa = \\frac{1}{2}\\nabla^2\\delta\\psi$.\n", + "\n", + "The visibility-space implementation uses the **sparse-operator (w-tilde) route**, whose cost scales with the\n", + "number of real-space mask pixels \u2014 independent of the number of visibilities \u2014 making it the path for real\n", + "datasets (ALMA / VLBI scale). This is the regime of the technique's benchmark detections: the ~$10^6 M_\\odot$\n", + "object in JVAS B1938+666 (Powell et al. 2025, Nat. Astron. 9, 1714; Vegetti et al. 2026).\n", + "\n", + "The workflow below is the configuration certified by the PyAutoLens validation campaign, which on this exact\n", + "setup recovers the simulated $10^{10} M_\\odot$ subhalo with dkappa correlation ~0.83 against the truth, the\n", + "peak ~0.15\" from the true position and ~6 sigma significance over the field:\n", + "\n", + " 1. verify the smooth-model source fit reaches $\\chi^2$/dof ~ 1-3 (the \"regime gate\" \u2014 corrections only\n", + " measure residual mass structure once the source model fits the data to near the noise);\n", + " 2. run the one-shot joint inversion (`FitDpsiSrcInterferometer`) with an arc-restricted `dpsi_mask`;\n", + " 3. refine with the iterative engine **warm-started from the one-shot**, with the regularization strengths\n", + " re-optimized by Bayesian evidence at every step (the Koopmans 2005 / Vegetti & Koopmans 2009 scheme).\n", + "\n", + "If you use this functionality in your research, please cite Cao et al. 2025, whose `potential_correction`\n", + "package (https://github.com/caoxiaoyue/lensing_potential_correction) the implementation is ported from;\n", + "citation materials are provided at https://github.com/caoxiaoyue/potential_correction_paper.\n", + "\n", + "This script takes a few minutes to run (the one-off sparse-operator precomputation dominates).\n", + "\n", + "__Prerequisites__\n", + "\n", + "- `guides/advanced/potential_correction.ipynb` \u2014 the technique overview (imaging).\n", + "- `interferometer/start_here.ipynb` \u2014 interferometer dataset basics." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Google Colab Setup__\n", + "\n", + "This cell sets up the environment when the notebook is run on Google Colab: it installs the\n", + "required PyAuto packages, clones the workspace (configuration files and example datasets) and\n", + "points the configuration at it. If you are running the notebook elsewhere (e.g. locally via\n", + "your own installation) it does nothing, and you can run it safely.\n", + "\n", + "Colab tip: model-fits run much faster on a GPU \u2014 enable one via \"Runtime\" -> \"Change runtime\n", + "type\" -> \"Hardware accelerator\" before running the notebook." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "try:\n", + " import google.colab\n", + " import subprocess\n", + " import sys\n", + "\n", + " subprocess.check_call(\n", + " [sys.executable, \"-m\", \"pip\", \"install\", \"autoconf\", \"--no-deps\"]\n", + " )\n", + "except ImportError:\n", + " pass\n", + "\n", + "from autoconf import setup_colab\n", + "\n", + "setup_colab.setup(\"autolens\")" + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "\n", + "# %matplotlib inline\n", + "# from pyprojroot import here\n", + "# workspace_path = str(here())\n", + "# %cd $workspace_path\n", + "# print(f\"Working Directory has been set to `{workspace_path}`\")\n", + "\n", + "import os\n", + "import numpy as np\n", + "\n", + "import autoarray as aa\n", + "import autolens as al" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__UV Coverage__\n", + "\n", + "We simulate a small but realistic dataset in-memory (seeded, so this script is fully reproducible). Random uv\n", + "points produce heavy sidelobes that smear the corrections, so we build an ALMA-like earth-rotation-synthesis\n", + "distribution: Gaussian-distributed antennas observed over an hour-angle track, giving the dense elliptical-ring\n", + "uv sampling of a real observation. The maximum baseline is kept at ~60% of the real-space grid's Nyquist limit.\n", + "\n", + "For real data, load your visibilities via `al.Interferometer.from_fits` as in `interferometer/start_here.py`." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "\n", + "\n", + "def synthesis_uv_from(n_ant, n_times, max_baseline_wavelengths, seed=0):\n", + " rng = np.random.default_rng(seed)\n", + " antennas = rng.normal(scale=max_baseline_wavelengths / 4.0, size=(n_ant, 2))\n", + " i_idx, j_idx = np.triu_indices(n_ant, k=1)\n", + " baselines = antennas[i_idx] - antennas[j_idx]\n", + "\n", + " hour_angles = np.linspace(-np.pi / 3.0, np.pi / 3.0, n_times)\n", + " declination = np.deg2rad(-23.0)\n", + "\n", + " uv_list = []\n", + " for ha in hour_angles:\n", + " u = baselines[:, 0] * np.cos(ha) - baselines[:, 1] * np.sin(ha)\n", + " v = (\n", + " baselines[:, 0] * np.sin(ha) * np.sin(declination)\n", + " + baselines[:, 1] * np.cos(ha) * np.sin(declination)\n", + " )\n", + " uv_list.append(np.stack([u, v], axis=1))\n", + " return np.concatenate(uv_list, axis=0)\n", + "\n", + "\n", + "pixel_scale = 0.08\n", + "nyquist = 1.0 / (2.0 * pixel_scale * np.pi / 180.0 / 3600.0)\n", + "uv_wavelengths = synthesis_uv_from(\n", + " n_ant=18, n_times=8, max_baseline_wavelengths=0.6 * nyquist\n", + ")\n", + "print(f\"{uv_wavelengths.shape[0]} visibilities\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Simulate__\n", + "\n", + "An `Isothermal` lens whose true mass also contains a $10^{10} M_\\odot$ NFW subhalo on the Einstein ring, and a\n", + "compact double-Gaussian source." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "real_space_mask = al.Mask2D.circular(\n", + " shape_native=(72, 72), pixel_scales=pixel_scale, radius=2.6\n", + ")\n", + "grid = al.Grid2D.from_mask(mask=real_space_mask)\n", + "\n", + "subhalo_centre = (1.41, 0.0)\n", + "true_subhalo = al.mp.NFWMCRLudlowSph(\n", + " centre=subhalo_centre, mass_at_200=1.0e10, redshift_object=0.2, redshift_source=0.6\n", + ")\n", + "lens_true = al.Galaxy(\n", + " redshift=0.2,\n", + " mass=al.mp.Isothermal(\n", + " centre=(0.0, 0.0),\n", + " einstein_radius=1.4,\n", + " ell_comps=al.convert.ell_comps_from(axis_ratio=0.9, angle=0.0),\n", + " ),\n", + " subhalo=true_subhalo,\n", + ")\n", + "source_true = al.Galaxy(\n", + " redshift=0.6,\n", + " bulge0=al.lp.Gaussian(\n", + " centre=(0.0, 0.0),\n", + " ell_comps=al.convert.ell_comps_from(axis_ratio=0.6, angle=45.0),\n", + " intensity=5.0,\n", + " sigma=0.15,\n", + " ),\n", + " bulge1=al.lp.Gaussian(\n", + " centre=(0.0, 0.4),\n", + " ell_comps=al.convert.ell_comps_from(axis_ratio=0.4, angle=135.0),\n", + " intensity=3.0,\n", + " sigma=0.1,\n", + " ),\n", + ")\n", + "\n", + "simulator = al.SimulatorInterferometer(\n", + " uv_wavelengths=uv_wavelengths,\n", + " exposure_time=300.0,\n", + " noise_sigma=2.0,\n", + " noise_seed=1,\n", + ")\n", + "dataset = simulator.via_tracer_from(\n", + " tracer=al.Tracer(galaxies=[lens_true, source_true]), grid=grid\n", + ")\n", + "dataset = al.Interferometer(\n", + " data=dataset.data,\n", + " noise_map=dataset.noise_map,\n", + " uv_wavelengths=uv_wavelengths,\n", + " real_space_mask=real_space_mask,\n", + ")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Sparse Operator__\n", + "\n", + "The w-tilde machinery is precomputed once per dataset. For large visibility counts this is the expensive step\n", + "(it is also noise-map dependent) \u2014 for repeated analyses of the same dataset, cache it to disk and reload via\n", + "`apply_sparse_operator(nufft_precision_operator=np.load(...))`." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "dataset = dataset.apply_sparse_operator()" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Smooth Starting Model__\n", + "\n", + "Potential corrections perturb a smooth starting model \u2014 in a real analysis, your maximum-likelihood parametric\n", + "fit (see `interferometer/modeling`). Here we use the true smooth lens (without the subhalo) and the true\n", + "source, isolating the technique itself. The source enters the correction operator through a `SrcFactory`: an\n", + "`AnalyticSrcFactory` here; use `PixSrcFactoryITP` to build it from a pixelized reconstruction of a previous fit." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "lens_smooth = al.Galaxy(redshift=0.2, mass=lens_true.mass)\n", + "source_start = al.pc.AnalyticSrcFactory(source_galaxy=source_true)" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Source Pixelization + Regime Gate__\n", + "\n", + "The source is reconstructed on a k-nearest-neighbor mesh distributed via an `Overlay` image-mesh. Before\n", + "trusting any correction, verify the smooth model + this source pixelization actually fit the data: if\n", + "$\\chi^2$/dof is far above ~1-3, the corrections will absorb source-model error instead of mass structure and\n", + "the dkappa map is meaningless (this gate is the single most common failure mode in practice)." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "grid_slim = dataset.grid.slim\n", + "source_shape = (\n", + " int(float(grid_slim[:, 0].max() - grid_slim[:, 0].min()) / pixel_scale / 2.0),\n", + " int(float(grid_slim[:, 1].max() - grid_slim[:, 1].min()) / pixel_scale / 2.0),\n", + ")\n", + "src_pixelization = al.Pixelization(\n", + " mesh=al.mesh.KNearestNeighbor(pixels=int(np.prod(source_shape))),\n", + " regularization=al.reg.Constant(coefficient=1.0),\n", + ")\n", + "src_image_mesh = al.image_mesh.Overlay(shape=source_shape)\n", + "\n", + "from autogalaxy.analysis.adapt_images.adapt_images import AdaptImages\n", + "\n", + "source_galaxy_pix = al.Galaxy(redshift=0.6, pixelization=src_pixelization)\n", + "image_plane_mesh_grid = src_image_mesh.image_plane_mesh_grid_from(mask=real_space_mask)\n", + "smooth_fit = al.FitInterferometer(\n", + " dataset=dataset,\n", + " tracer=al.Tracer(galaxies=[lens_smooth, source_galaxy_pix]),\n", + " adapt_images=AdaptImages(\n", + " galaxy_image_plane_mesh_grid_dict={source_galaxy_pix: image_plane_mesh_grid}\n", + " ),\n", + " settings=aa.Settings(use_positive_only_solver=True, use_border_relocator=True),\n", + ")\n", + "chi2_dof = float(smooth_fit.chi_squared) / (2 * uv_wavelengths.shape[0])\n", + "print(f\"regime gate: smooth-model chi2/dof = {chi2_dof:.2f} (target ~1-3)\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Arc-Restricted Dpsi Mesh__\n", + "\n", + "The real-space mask must stay a filled circle (it defines the sparse operator's FFT extent), but the\n", + "corrections are only constrained where the lensed arcs are \u2014 restrict the dpsi mesh to an arc-tracing sub-mask\n", + "built from the smooth-model image geometry." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "tracer_smooth = al.Tracer(galaxies=[lens_smooth, source_true])\n", + "arc_image = np.asarray(tracer_smooth.image_2d_from(grid=grid).native)\n", + "arc_mask = al.pc.util.arc_mask_from(\n", + " arc_image / (0.05 * arc_image.max()), threshold=3.0, ignore_size=10, ext_size=3\n", + ")\n", + "dpsi_mask = ~((~arc_mask) & (~np.asarray(real_space_mask)))\n", + "print(f\"dpsi mesh restricted to {int(np.count_nonzero(~dpsi_mask))} arc pixels\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__One-Shot Joint Inversion__\n", + "\n", + "The pixelized source and dpsi are solved in one linear system through the sparse route, with Matern\n", + "regularization on the corrections (the Cao et al. 2025 scheme, which recovers localised perturbers where\n", + "curvature penalties smear them out)." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "dpsi_pixelization = al.pc.DpsiPixelization(\n", + " mesh=al.pc.RegularDpsiMesh(factor=2),\n", + " regularization=al.reg.MaternKernel(coefficient=2000.0, scale=4.0, nu=2.5),\n", + ")\n", + "\n", + "fit = al.pc.FitDpsiSrcInterferometer(\n", + " dataset=dataset,\n", + " lens_start=lens_smooth,\n", + " source_start=source_start,\n", + " dpsi_pixelization=dpsi_pixelization,\n", + " src_pixelization=src_pixelization,\n", + " src_image_mesh=src_image_mesh,\n", + " dpsi_mask=dpsi_mask,\n", + " use_sparse_operator=True,\n", + ")\n", + "print(f\"one-shot joint log evidence = {fit.log_evidence:.4e}\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Dkappa Recovery Metrics__\n", + "\n", + "Because this is a simulation we can compare the recovered convergence correction to the true subhalo's\n", + "convergence: the correlation over the dpsi mesh, the peak's offset from the true position, and the peak's\n", + "significance over the field far from the subhalo." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "\n", + "\n", + "def dkappa_metrics(pair_obj, dkappa_rec, tag):\n", + " points = np.vstack([pair_obj.ygrid_dpsi_1d, pair_obj.xgrid_dpsi_1d]).T\n", + " dkappa_true = np.asarray(\n", + " true_subhalo.convergence_2d_from(grid=al.Grid2DIrregular(values=points))\n", + " )\n", + " corr = float(np.corrcoef(dkappa_rec, dkappa_true)[0, 1])\n", + " peak = points[int(np.argmax(dkappa_rec))]\n", + " dist = float(np.hypot(peak[0] - subhalo_centre[0], peak[1] - subhalo_centre[1]))\n", + " r = np.hypot(points[:, 0] - subhalo_centre[0], points[:, 1] - subhalo_centre[1])\n", + " significance = float(dkappa_rec.max() / dkappa_rec[r > 1.5].std())\n", + " print(\n", + " f\"{tag}: corr(dkappa_rec, dkappa_true) = {corr:.3f}; \"\n", + " f\"peak offset = {dist:.2f}\\\"; significance = {significance:.1f} sigma\"\n", + " )\n", + " return corr, dist\n", + "\n", + "\n", + "dkappa = np.asarray(fit.pair_dpsi_data_obj.hamiltonian_dpsi @ fit.best_fit_dpsi)\n", + "dkappa_metrics(fit.pair_dpsi_data_obj, dkappa, \"one-shot\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Iterative Refinement__\n", + "\n", + "The iterative Levenberg-Marquardt engine re-ray-traces through the corrected lens at each accepted step. Two\n", + "practices are essential (both certified by the validation campaign, mirroring the published pipelines):\n", + "\n", + "- **warm-start from the one-shot solution** (`x0=`), so the LM refines inside the right basin;\n", + "- **evidence-control the regularization strengths** (`reg_optimize_every=1`) rather than fixing them.\n", + "\n", + "Under JAX, pass `xp=jax.numpy` (with 64-bit enabled) to run the LM kernels JIT-compiled on GPU \u2014 the API is\n", + "identical, per the ecosystem's xp convention." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "iter_fit = al.pc.IterFitDpsiSrcInterferometer(\n", + " dataset=dataset,\n", + " lens_start=lens_smooth,\n", + " dpsi_pixelization=dpsi_pixelization,\n", + " src_pixelization=src_pixelization,\n", + " src_image_mesh=src_image_mesh,\n", + " dpsi_mask=dpsi_mask,\n", + " gauge_constraints=True,\n", + " n_iter=6,\n", + " reg_optimize_every=1,\n", + ")\n", + "s_opt, dpsi_opt = iter_fit.solve_joint_optimization(x0=np.asarray(fit.src_dpsi_slim))\n", + "print(f\"iterative Laplace log evidence = {iter_fit.log_evidence(s=s_opt, dpsi=dpsi_opt):.4e}\")\n", + "\n", + "dkappa_iter = np.asarray(iter_fit.pair_dpsi_data_obj.hamiltonian_dpsi @ dpsi_opt)\n", + "dkappa_metrics(iter_fit.pair_dpsi_data_obj, dkappa_iter, \"iterative\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Dkappa Map__\n", + "\n", + "The convergence-correction map \u2014 the subhalo appears as a positive peak at its position (1.41\", 0.0\")." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "os.makedirs(\"output\", exist_ok=True)\n", + "fig, axes = plt.subplots(1, 2, figsize=(10, 4))\n", + "al.pc.visualize.imshow_masked_data(\n", + " dkappa,\n", + " fit.pair_dpsi_data_obj.mask_dpsi,\n", + " ax=axes[0],\n", + " origin=\"upper\",\n", + " extent=fit.pair_dpsi_data_obj.data_bound,\n", + ")\n", + "axes[0].set_title(\"dkappa (one-shot)\")\n", + "al.pc.visualize.imshow_masked_data(\n", + " dkappa_iter,\n", + " iter_fit.pair_dpsi_data_obj.mask_dpsi,\n", + " ax=axes[1],\n", + " origin=\"upper\",\n", + " extent=iter_fit.pair_dpsi_data_obj.data_bound,\n", + ")\n", + "axes[1].set_title(\"dkappa (iterative)\")\n", + "plt.tight_layout()\n", + "plt.savefig(\n", + " os.path.join(\"output\", \"potential_correction_interferometer.png\"),\n", + " bbox_inches=\"tight\",\n", + ")\n", + "plt.close()" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Wrap Up__\n", + "\n", + "For hyper-parameter sampling with a non-linear search, wrap the fits in `al.pc.DpsiSrcInvInterferometerAnalysis`\n", + "(one-shot) or `al.pc.IterDpsiSrcInvInterferometerAnalysis` (iterative) \u2014 the Bayesian evidence consistently\n", + "ranks over-fit and degenerate solutions below recovering ones, so evidence-driven model selection self-protects.\n", + "\n", + "If you use this functionality, please cite Cao et al. 2025\n", + "(https://github.com/caoxiaoyue/potential_correction_paper) alongside PyAutoLens." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [], + "outputs": [], + "execution_count": null + } + ], + "metadata": { + "anaconda-cloud": {}, + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.1" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} \ No newline at end of file diff --git a/notebooks/weak/features/strong_lensing/a2744.ipynb b/notebooks/weak/features/strong_lensing/a2744.ipynb index 865fe2841..f8620984d 100644 --- a/notebooks/weak/features/strong_lensing/a2744.ipynb +++ b/notebooks/weak/features/strong_lensing/a2744.ipynb @@ -230,9 +230,11 @@ "galaxy_models = al.galaxy_af_models_from_csv_tables(mass_table, point_table)\n", "\n", "for name in (\"lens_0\", \"lens_1\"):\n", - " galaxy_models[name].mass.ra = af.UniformPrior(lower_limit=1.0, upper_limit=15.0)\n", - " galaxy_models[name].mass.rs = af.UniformPrior(lower_limit=5.0, upper_limit=40.0)\n", - " galaxy_models[name].mass.b0 = af.UniformPrior(lower_limit=0.1, upper_limit=10.0)\n", + " galaxy_models[name].mass.sigma = af.UniformPrior(lower_limit=50.0, upper_limit=600.0)\n", + " galaxy_models[name].mass.r_core = af.UniformPrior(lower_limit=1.0, upper_limit=15.0)\n", + " galaxy_models[name].mass.r_cut = af.UniformPrior(lower_limit=5.0, upper_limit=40.0)\n", + " galaxy_models[name].mass.H0 = 67.66\n", + " galaxy_models[name].mass.Om0 = 0.30966\n", "\n", "galaxy_models[\"host_halo\"].dark.mass_at_200 = af.LogUniformPrior(\n", " lower_limit=10**14.5, upper_limit=10**16.0\n", @@ -248,11 +250,13 @@ " mean=float(np.mean(positions[:, 1])), sigma=3.0\n", " )\n", "\n", - "scaling_b0_ref = af.UniformPrior(lower_limit=0.0, upper_limit=1.0)\n", - "scaling_exponent = 0.5\n", + "scaling_sigma_ref = af.UniformPrior(lower_limit=0.0, upper_limit=300.0)\n", + "scaling_sigma_exponent = 0.25\n", + "scaling_radius_exponent = 0.5\n", "reference_luminosity = 1.0\n", - "scaling_ra_ref_fixed = 0.158\n", - "scaling_rs_ref_fixed = 15.8\n", + "scaling_r_core_ref_fixed = 0.158\n", + "scaling_r_cut_ref_fixed = 15.8\n", + "source_redshift_max = max(float(d.redshift) for d in dataset_list)\n", "\n", "scaling_galaxies_list = []\n", "for centre, luminosity in zip(\n", @@ -262,9 +266,13 @@ "\n", " mass = af.Model(al.mp.dPIEMassSph)\n", " mass.centre = tuple(centre)\n", - " mass.ra = scaling_ra_ref_fixed * luminosity_ratio**scaling_exponent\n", - " mass.rs = scaling_rs_ref_fixed * luminosity_ratio**scaling_exponent\n", - " mass.b0 = scaling_b0_ref * luminosity_ratio**scaling_exponent\n", + " mass.sigma = scaling_sigma_ref * luminosity_ratio**scaling_sigma_exponent\n", + " mass.r_core = scaling_r_core_ref_fixed * luminosity_ratio**scaling_radius_exponent\n", + " mass.r_cut = scaling_r_cut_ref_fixed * luminosity_ratio**scaling_radius_exponent\n", + " mass.redshift_object = redshift_lens\n", + " mass.redshift_source = source_redshift_max\n", + " mass.H0 = 67.66\n", + " mass.Om0 = 0.30966\n", "\n", " scaling_galaxies_list.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass))\n", "\n", diff --git a/scripts/cluster/csv_api.py b/scripts/cluster/csv_api.py index 38451be2c..f76ba0e36 100644 --- a/scripts/cluster/csv_api.py +++ b/scripts/cluster/csv_api.py @@ -193,8 +193,26 @@ # family carries the 2 source-galaxy Point components. mass_profiles = { - "lens_0": {"mass": al.mp.dPIEMassSph(centre=(0.0, 0.0), ra=8.0, rs=20.0, b0=3.0)}, - "lens_1": {"mass": al.mp.dPIEMassSph(centre=(10.0, 8.0), ra=5.0, rs=12.0, b0=1.2)}, + "lens_0": { + "mass": al.mp.dPIEMassSph( + centre=(0.0, 0.0), + sigma=330.0, + r_core=8.0, + r_cut=20.0, + redshift_object=redshift_lens, + redshift_source=max(source_redshifts), + ) + }, + "lens_1": { + "mass": al.mp.dPIEMassSph( + centre=(10.0, 8.0), + sigma=210.0, + r_core=5.0, + r_cut=12.0, + redshift_object=redshift_lens, + redshift_source=max(source_redshifts), + ) + }, "host_halo": { "dark": al.mp.NFWMCRLudlowSph( centre=(0.0, 0.0), @@ -348,9 +366,9 @@ mutate the returned model: ```python -galaxy_models["lens_0"].mass.ra = af.UniformPrior(lower_limit=1.0, upper_limit=15.0) -galaxy_models["lens_0"].mass.rs = af.UniformPrior(lower_limit=5.0, upper_limit=40.0) -galaxy_models["lens_0"].mass.b0 = af.UniformPrior(lower_limit=0.1, upper_limit=10.0) +galaxy_models["lens_0"].mass.sigma = af.UniformPrior(lower_limit=50.0, upper_limit=600.0) +galaxy_models["lens_0"].mass.r_core = af.UniformPrior(lower_limit=1.0, upper_limit=15.0) +galaxy_models["lens_0"].mass.r_cut = af.UniformPrior(lower_limit=5.0, upper_limit=40.0) ``` This is the same composition pattern as ``af.Model(al.Galaxy, mass=...)`` @@ -363,11 +381,12 @@ point_table, ) -# Mutate selected params on the main-lens mass profiles into priors. +# Mutate selected params on the main-lens mass profiles into priors (sigma is +# Lenstool's fiducial v_disp in km/s; radii in arcsec). for name in ("lens_0", "lens_1"): - galaxy_models[name].mass.ra = af.UniformPrior(lower_limit=1.0, upper_limit=15.0) - galaxy_models[name].mass.rs = af.UniformPrior(lower_limit=5.0, upper_limit=40.0) - galaxy_models[name].mass.b0 = af.UniformPrior(lower_limit=0.1, upper_limit=10.0) + galaxy_models[name].mass.sigma = af.UniformPrior(lower_limit=50.0, upper_limit=600.0) + galaxy_models[name].mass.r_core = af.UniformPrior(lower_limit=1.0, upper_limit=15.0) + galaxy_models[name].mass.r_cut = af.UniformPrior(lower_limit=5.0, upper_limit=40.0) # Host halo: free mass_at_200, keep centre + redshifts fixed. galaxy_models["host_halo"].dark.mass_at_200 = af.LogUniformPrior( @@ -391,7 +410,7 @@ member and emitting an ``attr_name`` column would be more overhead than signal — every row uses the same ``dPIEMassSph`` mass profile with parameters derived from the reference-anchored scaling relation's shared -``b0_ref`` normalization (see ``modeling.py``). +``sigma_ref`` normalization (see ``modeling.py``). ``al.galaxy_table_to_csv`` and ``al.galaxy_table_from_csv`` are the schema-specific writers/readers. The simulator emits 10 scaling members @@ -439,18 +458,18 @@ __Lenstool-Parameterized Rows__ -Because the ``profile_class`` column dispatches against the full ``al.mp`` namespace, a -``mass.csv`` can carry rows in **Lenstool's native parameterization** via ``dPIEMassLenstool`` -— the columns become the ``.par``-file keywords verbatim:: +The default dPIE classes ARE Lenstool's native parameterization, so an elliptical ``dPIEMass`` +row's columns are the ``.par``-file keywords verbatim:: galaxy,attr_name,profile_class,y,x,ellipticity,angle_pos,sigma,r_core,r_cut,redshift_object,redshift_source,H0,Om0,redshift - O1,mass,dPIEMassLenstool,1.479,-2.997,0.678,8.971,987.34,18.96,283.54,0.39,11.76,70.0,0.3,0.39 + O1,mass,dPIEMass,1.479,-2.997,0.678,8.971,987.34,18.96,283.54,0.39,11.76,70.0,0.3,0.39 ``sigma`` is Lenstool's fiducial ``v_disp`` (sigma_LT), radii are in arcsec, and the run's own cosmology travels as the flat ``H0`` / ``Om0`` columns. ``scripts/cluster/lenstool/`` builds its entire 149-component published model this way — the ``.par`` file becomes one canonical CSV. Note the multi-plane convention: ``redshift_source`` must be the tracer's *final* (highest) source -plane. +plane. (The internal ``(ra, rs, b0)`` parameterization remains available for CSV rows via the +non-standard ``dPIEMassB0`` / ``dPIEMassB0Sph`` classes.) Light-profile CSVs (``light.csv``) support the linear / operated variants with qualified class names (``linear.Sersic``, ``operated.Gaussian``); plain names resolve to the standard profiles. diff --git a/scripts/cluster/lenstool/README.md b/scripts/cluster/lenstool/README.md index 6d5c5193c..cd1987c82 100644 --- a/scripts/cluster/lenstool/README.md +++ b/scripts/cluster/lenstool/README.md @@ -43,7 +43,7 @@ The worked example is **SMACS J0723.3−7327**, the first JWST cluster, using th (median RMS 0.07″, all below 0.29″ — consistent with the published image-plane RMS of 0.32″ through typical magnifications), and the optional Verification II forward-solve checks the image-plane RMS directly. One convention this exercise pinned down the honest way: PyAutoLens's - multi-plane tracer normalizes profile deflections to its **final** plane, so `from_lenstool` + multi-plane tracer normalizes profile deflections to its **final** plane, so `dPIEMass` must be given the highest source redshift in the system as `redshift_source`. ## What's different (and why you might care) diff --git a/scripts/cluster/lenstool/data.py b/scripts/cluster/lenstool/data.py index b62669f78..66c47b3c5 100644 --- a/scripts/cluster/lenstool/data.py +++ b/scripts/cluster/lenstool/data.py @@ -21,7 +21,7 @@ - ``galcat.cat`` → ``members.csv`` (member catalogue: centres + shape/mag properties, the ``al.galaxy_table_from_csv`` schema) - ``best.par`` → ``mass.csv`` (every optimized ``potential`` section as one - ``dPIEMassLenstool`` row of the canonical named-galaxy model CSV — **the .par file as a + ``dPIEMass`` row of the canonical named-galaxy model CSV — **the .par file as a table**, read back with ``al.galaxy_models_from_csv`` like every other cluster dataset) __Attribution__ @@ -145,7 +145,7 @@ def lenstool_yx_from(ra: float, dec: float) -> tuple: scaling relation) and 144 cluster members whose parameters Lenstool derived from the scaling relation. Every section carries the *same five numbers you would type into PyAutoLens*: - Lenstool ``potential`` → ``dPIEMass.from_lenstool`` argument + Lenstool ``potential`` → ``dPIEMass`` constructor argument --------------------------------------------------------------------- x_centre / y_centre [arcsec] → centre=(y, x) ellipticity (a²-b²)/(a²+b²) → ellipticity @@ -268,7 +268,7 @@ def parse_best_par(path: Path) -> tuple: ) # The whole optimized model — 5 named halos + 144 scaling members — becomes ONE canonical -# ``mass.csv``: each ``potential`` section is a ``dPIEMassLenstool`` row whose columns are the +# ``mass.csv``: each ``potential`` section is a ``dPIEMass`` row whose columns are the # ``.par`` keywords verbatim (sigma, r_core, r_cut, ellipticity, angle_pos) plus the run's # redshifts and cosmology as flat values. ``modeling.py`` reads it back with the same # ``al.galaxy_models_from_csv`` call used throughout ``scripts/cluster/``. @@ -281,7 +281,7 @@ def parse_best_par(path: Path) -> tuple: for h in halos: name = h["label"] if h["label"].startswith("O") else f"member_{h['label']}" profiles_by_galaxy[name] = { - "mass": al.mp.dPIEMassLenstool( + "mass": al.mp.dPIEMass( centre=(h["y"], h["x"]), ellipticity=h["ellipticity"], angle_pos=h["angle_pos"], @@ -302,7 +302,7 @@ def parse_best_par(path: Path) -> tuple: redshifts={name: 0.39 for name in profiles_by_galaxy}, ) -print("Wrote point_datasets.csv, members.csv, mass.csv (149 dPIEMassLenstool rows).") +print("Wrote point_datasets.csv, members.csv, mass.csv (149 dPIEMass rows).") """ __Image Cutout__ diff --git a/scripts/cluster/lenstool/modeling.py b/scripts/cluster/lenstool/modeling.py index 18b81bc1e..0193b1509 100644 --- a/scripts/cluster/lenstool/modeling.py +++ b/scripts/cluster/lenstool/modeling.py @@ -39,7 +39,7 @@ - **sigma is the fiducial velocity dispersion sigma_LT**, not the physical central velocity dispersion: sigma_0 = sqrt(3/2) * sigma_LT (Eliasdottir et al. 2007, App. A). PyAutoLens's - ``from_lenstool`` / ``dPIEMassLenstool`` take sigma_LT — quote ``v_disp`` values unchanged. + ``dPIEMass`` takes sigma_LT — quote ``v_disp`` values unchanged. Feeding a *measured* stellar velocity dispersion here overestimates the mass by 50%. - **The x axis points West** in Lenstool's relative frame (data.py verifies this against the data). All coordinates in this script live in that frame, so numbers compare directly to the @@ -57,7 +57,7 @@ __Contents__ - **Load Data:** the CSVs written by ``data.py``. -- **The Published Model, Reconstructed:** 149 ``from_lenstool`` profiles + 21 point sources. +- **The Published Model, Reconstructed:** 149 Lenstool-native ``dPIEMass`` profiles + 21 point sources. - **Verification I — source-plane compactness:** observed images trace to tight source groups. - **Verification II — image-plane RMS (optional):** forward-solve vs the published 0.32". - **Critical Curves (optional):** the per-source-plane critical curves over the HST image. @@ -105,7 +105,7 @@ exist, the model-optimized values of ``best.par`` otherwise). - ``mass.csv`` — the complete optimized mass model in the **canonical named-galaxy CSV** (the same ``al.galaxy_models_from_csv`` format every cluster script uses): 149 rows of - ``profile_class = dPIEMassLenstool``, one per ``potential`` section of ``best.par``, whose + ``profile_class = dPIEMass``, one per ``potential`` section of ``best.par``, whose columns are the ``.par`` keywords verbatim. The five individually-optimized halos are named O1 (cluster-scale), O2 (BCG), O3 ("dNW"), O4 ("ICL"), O5 ("eCM"); the 144 scaling members are ``member_``. @@ -124,14 +124,15 @@ members_table = al.galaxy_table_from_csv(file_path=dataset_path / "members.csv") print( - f"mass.csv: {len(mass_table.rows)} dPIEMassLenstool rows | " + f"mass.csv: {len(mass_table.rows)} dPIEMass rows | " f"members.csv: {len(members_table.luminosities)} catalogue members" ) """ __The Published Model, Reconstructed__ -Every ``potential`` section of ``best.par`` becomes one ``al.mp.dPIEMass`` via ``from_lenstool`` +Every ``potential`` section of ``best.par`` becomes one ``al.mp.dPIEMass`` — the default dPIE is +Lenstool's native parameterization, so the ``.par`` keywords are the constructor arguments — the arguments are the ``.par`` keywords, verbatim. This is the whole point of the Lenstool-native API: nothing is transcribed by hand, and the sqrt(3/2) sigma convention, the ellipticity conversion and the D_LS/D_S normalization are handled (and unit-tested) inside PyAutoLens. @@ -147,7 +148,7 @@ """ Z_REF_SOURCE = max(float(dataset.redshift) for dataset in dataset_list) -# One call: every mass.csv row instantiates its dPIEMassLenstool with the .par values — +# One call: every mass.csv row instantiates its dPIEMass with the .par values — # the redshift_source (final-plane) normalization and the run's H0/Om0 travel inside the # CSV columns, so nothing here needs to remember them. lens_galaxies = list(al.galaxies_from_csv_tables(mass_table).values()) @@ -283,7 +284,7 @@ Everything above used the published answer. A real analysis *fits*: the model below reproduces the composition Lenstool optimized (``input.par``), using the Lenstool-parameterized profile -``al.mp.dPIEMassLenstool`` so every free parameter, prior bound and posterior number is in +``al.mp.dPIEMass`` so every free parameter, prior bound and posterior number is in Lenstool units: - **O1, cluster halo**: centre U(-5,5)" (both axes), ellipticity U(0,0.8), angle U(-90,90), @@ -296,7 +297,7 @@ r_cut U(50,1000) kpc, sigma U(0,700). [7 free] - **O5 "eCM"**: centre fixed; ellipticity U(0,0.6), angle U(-90,90), r_core U(0,10) kpc, r_cut U(10,200) kpc, sigma U(0.1,300). [5 free] - - **potfile members**: every catalogue member gets a ``dPIEMassLenstool`` with centre, shape and + - **potfile members**: every catalogue member gets a ``dPIEMass`` with centre, shape and angle *fixed to the light* (``galcat.cat``) and its sigma / r_cut derived from two shared free parameters exactly as ``potfile`` defines — @@ -329,7 +330,7 @@ members_table.properties["ellipticity"], members_table.properties["angle_pos"], ): - mass = af.Model(al.mp.dPIEMassLenstool) + mass = af.Model(al.mp.dPIEMass) mass.centre = tuple(centre) mass.ellipticity = ellipticity mass.angle_pos = angle_pos diff --git a/scripts/cluster/likelihood_function.py b/scripts/cluster/likelihood_function.py index 044234b9c..1dd85a78f 100644 --- a/scripts/cluster/likelihood_function.py +++ b/scripts/cluster/likelihood_function.py @@ -146,22 +146,28 @@ # REFERENCE_LUMINOSITY is an explicit fixed constant (Lenstool's "mag0"), not the # sample max, and matches the simulator truth so members are reproduced exactly. scaling_galaxies = [] -SCALING_B0_REF_TRUTH = 0.190 -SCALING_EXPONENT = 0.5 -SCALING_RA_REF = 0.158 -SCALING_RS_REF = 15.8 +SCALING_SIGMA_REF_TRUTH = 85.0 +SCALING_SIGMA_EXPONENT = 0.25 +SCALING_RADIUS_EXPONENT = 0.5 +SCALING_R_CORE_REF = 0.158 +SCALING_R_CUT_REF = 15.8 REFERENCE_LUMINOSITY = 1.0 for centre, luminosity in zip( scaling_table.centres.in_list, scaling_table.luminosities ): luminosity_ratio = luminosity / REFERENCE_LUMINOSITY - ra = SCALING_RA_REF * luminosity_ratio**SCALING_EXPONENT - b0 = SCALING_B0_REF_TRUTH * luminosity_ratio**SCALING_EXPONENT - rs = SCALING_RS_REF * luminosity_ratio**SCALING_EXPONENT scaling_galaxies.append( al.Galaxy( redshift=redshift_lens, - mass=al.mp.dPIEMassSph(centre=tuple(centre), ra=ra, rs=rs, b0=b0), + mass=al.mp.dPIEMassSph( + centre=tuple(centre), + sigma=SCALING_SIGMA_REF_TRUTH + * luminosity_ratio**SCALING_SIGMA_EXPONENT, + r_core=SCALING_R_CORE_REF * luminosity_ratio**SCALING_RADIUS_EXPONENT, + r_cut=SCALING_R_CUT_REF * luminosity_ratio**SCALING_RADIUS_EXPONENT, + redshift_object=redshift_lens, + redshift_source=max(source_redshifts), + ), ) ) @@ -171,8 +177,9 @@ The tracer carries: - 2 main lens galaxies (BCG + satellite) — individually-modelled dPIE mass profiles. - - 10 scaling-tier member galaxies — dPIE mass profiles whose ``ra``, ``rs`` and ``b0`` derive from the - reference-anchored scaling relation ``b0 = b0_ref × (L/L_ref)^0.5`` (Lenstool convention). + - 10 scaling-tier member galaxies — dPIE mass profiles whose ``sigma``, ``r_core`` and ``r_cut`` derive + from the reference-anchored scaling relation ``sigma = sigma_ref × (L/L_ref)^0.25``, + ``radii ∝ (L/L_ref)^0.5`` (Lenstool convention). - 1 host dark matter halo — ``NFWMCRLudlowSph`` at the cluster centre. - 2 source galaxies — ``Point`` profiles at distinct redshifts (multi-plane). diff --git a/scripts/cluster/modeling.py b/scripts/cluster/modeling.py index 02a6ae4d9..351dc2b40 100644 --- a/scripts/cluster/modeling.py +++ b/scripts/cluster/modeling.py @@ -43,8 +43,8 @@ - There are 2 main lens galaxies with ``dPIEMassSph`` total mass distributions, each with their centre fixed to the values written out by the simulator [6 parameters]. - - There are 10 scaling-tier member galaxies. Each carries a ``dPIEMassSph`` mass with centre fixed, - ``ra`` fixed, and ``b0`` / ``rs`` derived from the shared reference-anchored scaling relation and the + - There are 10 scaling-tier member galaxies. Each carries a ``dPIEMassSph`` mass with centre fixed and + ``sigma`` / ``r_core`` / ``r_cut`` derived from the shared reference-anchored scaling relation and the per-member luminosity [1 parameter total for the entire tier]. - There is 1 standalone ``NFWMCRLudlowSph`` host dark matter halo with its centre fixed and a free ``mass_at_200`` [1 parameter]. @@ -53,11 +53,16 @@ The number of free parameters and therefore the dimensionality of non-linear parameter space is N=12. +The ``dPIEMassSph`` profile is parameterized in Lenstool's native convention — free parameters +``sigma`` (fiducial velocity dispersion ``v_disp`` in km/s), ``r_core`` and ``r_cut`` (arcsec) — so the +fitted posterior reads like a Lenstool results table (see ``cluster/simulator.py`` __dPIE Mass Profile__ +for the conventions, and ``cluster/lenstool/`` for a published-model worked example). + The defining feature of cluster modeling is the scaling tier: 10 lower-mass members are fit jointly -with a *single free parameter* (``b0_ref``, the lens strength of a galaxy at the reference magnitude; -the relation's exponent is fixed at the Faber-Jackson value of 0.5, following the Lenstool convention). -Adding more members to ``scaling_galaxies.csv`` in the future does not grow the dimensionality of -parameter space. +with a *single free parameter* (``sigma_ref``, the fiducial velocity dispersion of a galaxy at the +reference magnitude; the relation's exponents are fixed at the Faber-Jackson values, following the +Lenstool convention). Adding more members to ``scaling_galaxies.csv`` in the future does not grow the +dimensionality of parameter space. __Simulation__ @@ -245,10 +250,10 @@ with a ``dPIEMassSph`` total mass profile whose centre is fixed to ``main_lens_centres[i]``. - ``scaling_galaxies``: The 10 scaling-tier cluster members. Each carries a ``dPIEMassSph`` mass with - centre fixed (from the CSV), ``ra`` fixed, and ``b0`` / ``rs`` derived from the shared - reference-anchored scaling relation (single free normalization ``b0_ref``, exponents fixed at 0.5) - plus the per-member luminosity. The whole tier contributes 1 free parameter to the model regardless - of how many members are in the CSV. + centre fixed (from the CSV) and ``sigma`` / ``r_core`` / ``r_cut`` derived from the shared + reference-anchored scaling relation (single free normalization ``sigma_ref``; exponents fixed at + 0.25 for sigma and 0.5 for the radii) plus the per-member luminosity. The whole tier contributes + 1 free parameter to the model regardless of how many members are in the CSV. - ``host_halo``: A single standalone ``Galaxy`` carrying the cluster's ``NFWMCRLudlowSph`` dark matter halo. The halo is *not* tied to any individual member — it sits "on top of" the members and @@ -282,12 +287,12 @@ We compose a lens model where: - The 2 main lens galaxies each have a ``dPIEMassSph`` mass profile with centre fixed and free - ``ra``, ``rs``, ``b0`` — 3 free parameters per galaxy [6 parameters]. - - The 10 scaling-tier members share a single free parameter: ``b0_ref``, the lens strength of a - galaxy at the reference magnitude. Each member's ``ra``, ``rs`` and ``b0`` are computed as - ``ra_ref * (L / L_ref) ** 0.5``, ``rs_ref * (L / L_ref) ** 0.5`` and ``b0_ref * (L / L_ref) ** 0.5`` - with the exponent fixed; ``ra_ref`` (0.158") and ``rs_ref`` (15.8") are held fixed at the simulator - truth values [1 parameter]. + ``sigma``, ``r_core``, ``r_cut`` — 3 free parameters per galaxy [6 parameters]. + - The 10 scaling-tier members share a single free parameter: ``sigma_ref``, the fiducial velocity + dispersion of a galaxy at the reference magnitude. Each member's parameters are computed as + ``sigma_ref * (L / L_ref) ** 0.25``, ``r_core_ref * (L / L_ref) ** 0.5`` and + ``r_cut_ref * (L / L_ref) ** 0.5`` with the exponents fixed; ``r_core_ref`` (0.158") and + ``r_cut_ref`` (15.8") are held fixed at the simulator truth values [1 parameter]. - The host halo has an ``NFWMCRLudlowSph`` mass profile with centre fixed and a free ``mass_at_200`` [1 parameter]. - Each source has a ``Point`` model with free ``centre_0`` / ``centre_1`` priors initialised from the @@ -301,30 +306,31 @@ published cluster strong-lensing analysis (Limousin et al. 2005; Eliasdottir et al. 2007; Bergamini et al. 2019): - ra_i = ra_ref * (L_i / L_ref) ** 0.5 - rs_i = rs_ref * (L_i / L_ref) ** 0.5 - b0_i = b0_ref * (L_i / L_ref) ** 0.5 + sigma_i = sigma_ref * (L_i / L_ref) ** 0.25 + r_core_i = r_core_ref * (L_i / L_ref) ** 0.5 + r_cut_i = r_cut_ref * (L_i / L_ref) ** 0.5 The reference luminosity ``L_ref`` is an **explicit fixed constant** (Lenstool's reference magnitude ``mag0``), *not* the maximum luminosity of the current sample. Anchoring to a fixed reference makes the normalization invariant to which galaxies are placed in the scaling tier, and gives the single free -parameter ``b0_ref`` a stable, interpretable meaning: the lens strength (roughly the Einstein radius) of +parameter ``sigma_ref`` a stable, interpretable meaning: the fiducial velocity dispersion (in km/s) of a galaxy *at the reference magnitude*, for which a prior range is easy to motivate — unlike an abstract multiplicative factor whose units depend on the (arbitrary) luminosity normalization. In a real analysis set ``L_ref`` to the BCG magnitude (or a characteristic L*); here we use a fiducial ``L_ref = 1.0``. Only luminosity *ratios* enter, so the CSV's luminosity units are irrelevant; magnitude catalogues convert via ``L_i / L_ref = 10 ** (0.4 * (m_ref - m_i))``. -The exponent is *fixed* at 0.5 rather than fitted: b0 ∝ sigma² for the dPIE, and Faber-Jackson -(sigma ∝ L^(1/4)) gives b0 ∝ L^(1/2). The core radius ``ra`` and truncation radius ``rs`` scale with the -same fixed exponent, mirroring Lenstool's r_core ∝ L^(1/2) and r_cut ∝ L^(1/2). Freeing the exponent (or -``rs_ref`` / ``ra_ref``) is a one-line change shown in the code comment below — useful as a systematics -test, at the cost of the degeneracy between normalization and slope that the fixed-exponent convention -exists to avoid. When member velocity dispersions are available, the standard refinement is to calibrate -the exponents kinematically (Bergamini et al. 2019: sigma ∝ L^0.27-0.28 from MUSE member kinematics, i.e. -a b0 exponent ≈ 0.55, with the r_cut exponent from the fundamental plane). - -The simulator's truth value is ``b0_ref = 0.190`` arcsec (at ``L_ref = 1.0``). The prior below is much +The exponents are *fixed* rather than fitted: Faber-Jackson (L ∝ sigma^4) gives sigma ∝ L^(1/4), and +constant mass-to-light (M ∝ sigma^2 * r_cut ∝ L) gives r_cut ∝ L^(1/2) — exactly Lenstool's ``potfile`` +scalings (``vdslope 4``, ``slope 4``), applied to the core radius too. Since the dPIE lens strength obeys +b0 ∝ sigma^2, this is equivalent to the b0 ∝ L^(1/2) scaling of the internal parameterization. Freeing an +exponent (or ``r_cut_ref`` / ``r_core_ref``) is a one-line change shown in the code comment below — +useful as a systematics test, at the cost of the degeneracy between normalization and slope that the +fixed-exponent convention exists to avoid. When member velocity dispersions are available, the standard +refinement is to calibrate the exponents kinematically (Bergamini et al. 2019: sigma ∝ L^0.27-0.28 from +MUSE member kinematics, with the r_cut exponent from the fundamental plane). + +The simulator's truth value is ``sigma_ref = 85.0`` km/s (at ``L_ref = 1.0``). The prior below is much wider than the truth to give the search room. """ redshift_lens = 0.5 @@ -337,11 +343,16 @@ galaxy_models = al.galaxy_af_models_from_csv_tables(mass_table, point_table) -# Main Lens Galaxies: free dPIE ra / rs / b0 on each; centre stays fixed at the CSV value. +# Main Lens Galaxies: free dPIE sigma / r_core / r_cut on each; centre and redshifts +# stay fixed at the CSV values, and the cosmology constants H0 / Om0 are pinned (they +# are model *constants*, not parameters to sample — if left unset they would inherit +# the config's default priors and float). for name in ("lens_0", "lens_1"): - galaxy_models[name].mass.ra = af.UniformPrior(lower_limit=1.0, upper_limit=15.0) - galaxy_models[name].mass.rs = af.UniformPrior(lower_limit=5.0, upper_limit=40.0) - galaxy_models[name].mass.b0 = af.UniformPrior(lower_limit=0.1, upper_limit=10.0) + galaxy_models[name].mass.sigma = af.UniformPrior(lower_limit=50.0, upper_limit=600.0) + galaxy_models[name].mass.r_core = af.UniformPrior(lower_limit=1.0, upper_limit=15.0) + galaxy_models[name].mass.r_cut = af.UniformPrior(lower_limit=5.0, upper_limit=40.0) + galaxy_models[name].mass.H0 = 67.66 + galaxy_models[name].mass.Om0 = 0.30966 # Host Halo: free mass_at_200; centre + redshift_object + redshift_source stay fixed. galaxy_models["host_halo"].dark.mass_at_200 = af.LogUniformPrior( @@ -363,7 +374,7 @@ mean=float(np.mean(positions[:, 1])), sigma=3.0 ) -# Scaling Tier Members (dPIEMassSph; ra, rs and b0 all derived from the +# Scaling Tier Members (dPIEMassSph; sigma, r_core and r_cut all derived from the # reference-anchored scaling relation). # # The reference luminosity is an EXPLICIT FIXED constant (Lenstool's reference @@ -372,23 +383,25 @@ # magnitude in a real analysis; here it is a fiducial L* = 1.0. reference_luminosity = 1.0 -# b0_ref is defined ONCE outside the loop — the tier's only free parameter, the -# lens strength of a reference-magnitude galaxy. Every member's ra, rs and b0 are -# derived by scaling the fixed reference values (ra_ref, rs_ref) or the free b0_ref -# by its luminosity ratio, with the exponent fixed at the Faber-Jackson value 0.5. -# The entire tier therefore contributes 1 free parameter regardless of how many +# sigma_ref is defined ONCE outside the loop — the tier's only free parameter, the +# fiducial velocity dispersion of a reference-magnitude galaxy (km/s, so the prior +# range below is physically interpretable). Every member's sigma, r_core and r_cut +# are derived by scaling the free sigma_ref or the fixed reference radii by its +# luminosity ratio, with the exponents fixed at the Faber-Jackson values +# (sigma ∝ L^0.25, radii ∝ L^0.5 — Lenstool's potfile vdslope/slope = 4). The +# entire tier therefore contributes 1 free parameter regardless of how many # members are in scaling_galaxies.csv. # -# To free the exponent as a systematics test, replace the fixed value with e.g. -# `scaling_exponent = af.UniformPrior(lower_limit=0.0, upper_limit=1.0)` — every -# member's b0 then derives from two shared parameters, as in older versions of -# this example. +# To free an exponent as a systematics test, replace the fixed value with e.g. +# `scaling_sigma_exponent = af.UniformPrior(lower_limit=0.0, upper_limit=0.5)` — +# every member's sigma then derives from two shared parameters. -scaling_b0_ref = af.UniformPrior(lower_limit=0.0, upper_limit=1.0) -scaling_exponent = 0.5 +scaling_sigma_ref = af.UniformPrior(lower_limit=0.0, upper_limit=200.0) +scaling_sigma_exponent = 0.25 +scaling_radius_exponent = 0.5 -scaling_ra_ref_fixed = 0.158 -scaling_rs_ref_fixed = 15.8 +scaling_r_core_ref_fixed = 0.158 +scaling_r_cut_ref_fixed = 15.8 scaling_galaxies_list = [] for centre, luminosity in zip( @@ -398,9 +411,13 @@ mass = af.Model(al.mp.dPIEMassSph) mass.centre = tuple(centre) - mass.ra = scaling_ra_ref_fixed * luminosity_ratio**scaling_exponent - mass.rs = scaling_rs_ref_fixed * luminosity_ratio**scaling_exponent - mass.b0 = scaling_b0_ref * luminosity_ratio**scaling_exponent + mass.sigma = scaling_sigma_ref * luminosity_ratio**scaling_sigma_exponent + mass.r_core = scaling_r_core_ref_fixed * luminosity_ratio**scaling_radius_exponent + mass.r_cut = scaling_r_cut_ref_fixed * luminosity_ratio**scaling_radius_exponent + mass.redshift_object = redshift_lens + mass.redshift_source = max(source_redshifts) + mass.H0 = 67.66 + mass.Om0 = 0.30966 scaling_galaxies_list.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass)) diff --git a/scripts/cluster/simulator.py b/scripts/cluster/simulator.py index 928d68fb8..f57cfacc8 100644 --- a/scripts/cluster/simulator.py +++ b/scripts/cluster/simulator.py @@ -67,9 +67,9 @@ relation. Each member is individually weak compared to the main galaxies or the host halo, but the population together perturbs the deflection field non-trivially — exactly the regime in which the scaling-relation tier of the modeling API earns its keep. The number of free parameters does not grow - with the number of scaling members; a single shared normalization `b0_ref` (the lens strength of a - galaxy at the reference magnitude, with the relation's exponent fixed at the Faber-Jackson value) - determines every member's mass from its luminosity. + with the number of scaling members; a single shared normalization `sigma_ref` (the fiducial velocity + dispersion of a galaxy at the reference magnitude, with the relation's exponents fixed at the + Faber-Jackson values) determines every member's mass from its luminosity. - `host_halo_galaxy`: A standalone `Galaxy` holding the cluster's `NFWMCRLudlowSph` dark matter halo. It is not tied to any individual member galaxy — the halo is a separate mass component sitting "on top of" @@ -87,49 +87,62 @@ The cluster member galaxies use the dual Pseudo-Isothermal Elliptical (dPIE) mass profile introduced in Eliasdottir 2007 (https://arxiv.org/abs/0710.5636), the de facto standard for cluster strong lens modeling. -In spherical form (`dPIEMassSph`), its parameters are: - - - `ra` (arcsec): the core radius, below which the density profile flattens (kept small, ~0.05–0.1" at z=0.5). - - `rs` (arcsec): the truncation radius, above which the density falls as R^-4 (kept ~10–30" for cluster members). - - `b0` (arcsec): the mass normalization, roughly setting the galaxy-scale Einstein radius. +In spherical form (`dPIEMassSph`), its parameters are **Lenstool's native ones** — the same numbers that +appear in a Lenstool ``.par`` file or a published cluster paper's results table: + + - `sigma` (km/s): the fiducial velocity dispersion (Lenstool's ``v_disp``, sigma_LT). This is *not* the + physical central dispersion — sigma_0 = sqrt(3/2) * sigma_LT (Eliasdottir et al. 2007, App. A) — so + quote ``.par``/table values unchanged and never feed a measured stellar dispersion here. + - `r_core` (arcsec): the core radius, below which the density profile flattens. + - `r_cut` (arcsec): the truncation radius, above which the density falls as R^-4 (~10–30" for cluster members). + - `redshift_object` / `redshift_source`: the redshifts entering the D_LS / D_S normalization of the + internal lens strength ``b0 = 6 * 648000 * (sigma / c)^2 * (D_LS / D_S)`` arcsec. In this multi-plane + cluster, ``redshift_source`` is anchored to the *furthest* source (``z = 2.0``) — the same convention + as the ``NFWMCRLudlowSph`` host halo below, because the tracer normalizes profile deflections to the + final source plane. + +The internal ``(ra, rs, b0)`` parameterization remains available as ``dPIEMassB0Sph`` (non-standard; b0 +is the lens strength in arcsec, with b0 proportional to sigma^2). Per-galaxy values for the 2 main-tier galaxies are hand-tuned below; for the 10 scaling-tier members they are derived from each member's luminosity via the relation described next. __Luminosity-Mass Scaling Relation__ -The 10 scaling members share a reference-anchored relation for the dPIE mass normalization — the -convention used by Lenstool and essentially every published cluster strong-lensing analysis -(Limousin et al. 2005; Eliasdottir et al. 2007; Bergamini et al. 2019): +The 10 scaling members share a reference-anchored relation — the convention used by Lenstool and +essentially every published cluster strong-lensing analysis (Limousin et al. 2005; Eliasdottir et al. +2007; Bergamini et al. 2019): - ra_i = ra_ref * (L_i / L_ref) ** 0.5 - rs_i = rs_ref * (L_i / L_ref) ** 0.5 - b0_i = b0_ref * (L_i / L_ref) ** 0.5 + sigma_i = sigma_ref * (L_i / L_ref) ** 0.25 + r_core_i = r_core_ref * (L_i / L_ref) ** 0.5 + r_cut_i = r_cut_ref * (L_i / L_ref) ** 0.5 where ``L_ref`` is an **explicit fixed reference luminosity** (Lenstool's reference magnitude ``mag0``), -*not* the maximum luminosity of the sample, and ``b0_ref`` is the lens strength of a galaxy *at that -reference magnitude*. Anchoring to a fixed reference makes the normalization physically interpretable — -the Einstein-radius-like strength of a reference-magnitude galaxy — and invariant to which galaxies are -placed in the tier, which is what makes a sensible prior range easy to define. The exponent is **fixed -at 0.5** rather than fitted: for the dPIE, ``b0`` is proportional to the velocity dispersion squared, and -the Faber-Jackson relation (L ∝ sigma^4, i.e. sigma ∝ L^(1/4)) then gives b0 ∝ L^(1/2). Lenstool applies -the same fixed-exponent scaling to the core and truncation radii (r_core ∝ L^(1/2), r_cut ∝ L^(1/2)), -which is why ``ra`` and ``rs`` scale here too. - -Truth values used in this simulator are ``b0_ref = 0.190`` arcsec, ``rs_ref = 15.8`` arcsec and -``ra_ref = 0.158`` arcsec at a fiducial reference luminosity ``L_ref = 1.0`` (the model in -``modeling.py`` / ``start_here.py`` uses the SAME reference, so members stay consistent by construction). Luminosities are log-spaced across roughly 0.05–0.40, so -per-member ``b0`` values run from ~0.042 to 0.12 arcsec — each member is individually well below the -BCG (``b0 = 3.0``) but the 10 of them together perturb the deflection field by ~10–15%. - -The modeling script promotes ``b0_ref`` to the tier's single free parameter and recovers the truth value -when fit to the simulated point datasets. Adding more scaling members amounts to adding rows to +*not* the maximum luminosity of the sample, and ``sigma_ref`` is the fiducial velocity dispersion of a +galaxy *at that reference magnitude*. Anchoring to a fixed reference makes the normalization physically +interpretable — the velocity dispersion of a reference-magnitude galaxy, in km/s — and invariant to which +galaxies are placed in the tier, which is what makes a sensible prior range easy to define. The exponents +are **fixed** rather than fitted: the Faber-Jackson relation (L ∝ sigma^4) gives sigma ∝ L^(1/4), and +constant mass-to-light (M ∝ sigma^2 * r_cut ∝ L) then gives r_cut ∝ L^(1/2) — Lenstool's ``potfile`` +applies exactly these scalings (``vdslope 4``, ``slope 4``), including to the core radius. Because the +lens strength b0 is proportional to sigma^2, this is equivalent to the b0 ∝ L^(1/2) scaling of the +internal parameterization. + +Truth values used in this simulator are ``sigma_ref = 85.0`` km/s, ``r_cut_ref = 15.8`` arcsec and +``r_core_ref = 0.158`` arcsec at a fiducial reference luminosity ``L_ref = 1.0`` (the model in +``modeling.py`` / ``start_here.py`` uses the SAME reference, so members stay consistent by construction). +Luminosities are log-spaced across roughly 0.05–0.40, so per-member dispersions run from ~40 to ~68 km/s +— each member is individually well below the BCG (``sigma = 330`` km/s) but the 10 of them together +perturb the deflection field by ~10–15%. + +The modeling script promotes ``sigma_ref`` to the tier's single free parameter and recovers the truth +value when fit to the simulated point datasets. Adding more scaling members amounts to adding rows to ``scaling_galaxies.csv`` — the tier's free-parameter count stays at 1. Note that only the luminosity *ratios* ``L_i / L_ref`` enter the relation, so the units of the luminosity column are irrelevant; observational catalogues quoting magnitudes convert via ``L_i / L_ref = 10 ** (0.4 * (m_ref - m_i))``. Kinematic calibrations of the exponent exist for when higher fidelity is needed — Bergamini et al. 2019 -measure sigma ∝ L^0.27-0.28 from MUSE member kinematics (b0 exponent ≈ 0.55) and derive the r_cut -exponent from the fundamental plane — but 0.5 is the standard default. +measure sigma ∝ L^0.27-0.28 from MUSE member kinematics and derive the r_cut exponent from the +fundamental plane (beta_cut = gamma - 2*alpha + 1) — but 0.25 / 0.5 are the standard defaults. __NFWMCRLudlow Host Halo__ @@ -296,9 +309,9 @@ parameters representative of cluster members: a larger central BCG and one smaller satellite galaxy. """ main_lens_dpie_params = [ - # (ra, rs, b0) per galaxy — arcsec - (8.0, 20.0, 3.0), # BCG — strongest - (5.0, 12.0, 1.2), # satellite + # (r_core, r_cut, sigma) per galaxy — arcsec, arcsec, km/s + (8.0, 20.0, 330.0), # BCG — strongest + (5.0, 12.0, 210.0), # satellite ] main_lens_sersic_params = [ @@ -308,7 +321,7 @@ ] main_lens_galaxies = [] -for centre, (ra, rs, b0), (intensity, effective_radius, sersic_index) in zip( +for centre, (r_core, r_cut, sigma), (intensity, effective_radius, sersic_index) in zip( main_lens_centres, main_lens_dpie_params, main_lens_sersic_params ): bulge = al.lp.SersicSph( @@ -317,7 +330,14 @@ effective_radius=effective_radius, sersic_index=sersic_index, ) - mass = al.mp.dPIEMassSph(centre=centre, ra=ra, rs=rs, b0=b0) + mass = al.mp.dPIEMassSph( + centre=centre, + sigma=sigma, + r_core=r_core, + r_cut=r_cut, + redshift_object=redshift_lens, + redshift_source=max(source_redshifts), + ) main_lens_galaxies.append(al.Galaxy(redshift=redshift_lens, bulge=bulge, mass=mass)) """ @@ -325,19 +345,20 @@ The 10 cluster members modelled collectively via the luminosity-mass scaling relation (see the ``__Luminosity-Mass Scaling Relation__`` section of the module docstring). The simulator hardcodes the -truth value of ``b0_ref`` (the lens strength of a galaxy at the reference magnitude) and derives each -member's ``ra``, ``rs`` and ``b0`` from its luminosity ratio to the reference, with the exponent fixed -at the Faber-Jackson value of 0.5. The reference luminosity is an EXPLICIT FIXED constant (Lenstool's -``mag0``), NOT the sample max — here a fiducial ``L_ref = 1.0``; the model in ``modeling.py`` / -``start_here.py`` uses the SAME reference so members stay consistent by construction. Light profiles use -the per-member luminosity as the central intensity so the rendered image visibly traces the scaling-tier -population. -""" -scaling_b0_ref_truth = 0.190 -scaling_exponent = 0.5 +truth value of ``sigma_ref`` (the fiducial velocity dispersion of a galaxy at the reference magnitude) +and derives each member's ``sigma``, ``r_core`` and ``r_cut`` from its luminosity ratio to the reference, +with the exponents fixed at the Faber-Jackson values (sigma ∝ L^0.25; radii ∝ L^0.5). The reference +luminosity is an EXPLICIT FIXED constant (Lenstool's ``mag0``), NOT the sample max — here a fiducial +``L_ref = 1.0``; the model in ``modeling.py`` / ``start_here.py`` uses the SAME reference so members stay +consistent by construction. Light profiles use the per-member luminosity as the central intensity so the +rendered image visibly traces the scaling-tier population. +""" +scaling_sigma_ref_truth = 85.0 +scaling_sigma_exponent = 0.25 +scaling_radius_exponent = 0.5 reference_luminosity = 1.0 -scaling_ra_ref = 0.158 -scaling_rs_ref = 15.8 +scaling_r_core_ref = 0.158 +scaling_r_cut_ref = 15.8 scaling_galaxies = [] for centre, luminosity in zip(scaling_galaxies_centres, scaling_galaxies_luminosities): @@ -348,10 +369,14 @@ sersic_index=3.0, ) luminosity_ratio = luminosity / reference_luminosity - ra = scaling_ra_ref * luminosity_ratio**scaling_exponent - b0 = scaling_b0_ref_truth * luminosity_ratio**scaling_exponent - rs = scaling_rs_ref * luminosity_ratio**scaling_exponent - mass = al.mp.dPIEMassSph(centre=centre, ra=ra, rs=rs, b0=b0) + mass = al.mp.dPIEMassSph( + centre=centre, + sigma=scaling_sigma_ref_truth * luminosity_ratio**scaling_sigma_exponent, + r_core=scaling_r_core_ref * luminosity_ratio**scaling_radius_exponent, + r_cut=scaling_r_cut_ref * luminosity_ratio**scaling_radius_exponent, + redshift_object=redshift_lens, + redshift_source=max(source_redshifts), + ) scaling_galaxies.append(al.Galaxy(redshift=redshift_lens, bulge=bulge, mass=mass)) """ @@ -538,8 +563,8 @@ def jitted_solve(tracer, source_plane_coordinate): Scaling up a real cluster to a larger member population is then a CSV-level edit: add a row per additional member, fill in its centre and luminosity, save. The modeling script picks up the new rows -automatically and the scaling tier's free-parameter count stays at 1 (``b0_ref``, the reference -member's lens strength; the relation's exponents stay fixed at 0.5). Only luminosity *ratios* enter the +automatically and the scaling tier's free-parameter count stays at 1 (``sigma_ref``, the reference +member's fiducial velocity dispersion; the relation's exponents stay fixed). Only luminosity *ratios* enter the relation, so any consistent luminosity convention works — including converting from magnitudes. """ al.galaxy_table_to_csv( diff --git a/scripts/cluster/start_here.py b/scripts/cluster/start_here.py index a62284cb8..d1c3736ee 100644 --- a/scripts/cluster/start_here.py +++ b/scripts/cluster/start_here.py @@ -245,16 +245,18 @@ - **Main lens galaxies (2):** the two brightest core members (the BCG region galaxies), individually-modelled ``dPIEMassSph`` profiles with centre fixed to the observed light centres - and free ``ra``, ``rs``, ``b0``. **6 free parameters total.** - - - **Scaling-tier members (188):** ``dPIEMassSph`` profiles with centre fixed to the CSV centres. ``ra``, - ``rs`` and ``b0`` all derive from the reference-anchored relation used by Lenstool and standard in - published cluster analyses: ``ra = ra_ref * (L / L_ref) ** 0.5``, ``rs = rs_ref * (L / L_ref) ** 0.5`` - and ``b0 = b0_ref * (L / L_ref) ** 0.5``, where ``L_ref`` is an explicit fixed reference luminosity - (Lenstool's ``mag0``), *not* the sample max. The exponent is fixed at the Faber-Jackson value - (b0 ∝ sigma² and sigma ∝ L^(1/4) give b0 ∝ L^(1/2)) — only the normalization ``b0_ref``, the lens - strength of a reference-magnitude galaxy, is fitted. Our member luminosities are normalised to the - BCG's F160W flux, so ``L_ref = 1.0`` anchors the relation to the BCG itself. + and free ``sigma``, ``r_core``, ``r_cut`` — Lenstool's native dPIE parameters, so the posterior + reads like a Lenstool results table. **6 free parameters total.** + + - **Scaling-tier members (188):** ``dPIEMassSph`` profiles with centre fixed to the CSV centres. + ``sigma``, ``r_core`` and ``r_cut`` all derive from the reference-anchored relation used by Lenstool + and standard in published cluster analyses: ``sigma = sigma_ref * (L / L_ref) ** 0.25``, + ``r_core = r_core_ref * (L / L_ref) ** 0.5`` and ``r_cut = r_cut_ref * (L / L_ref) ** 0.5``, where + ``L_ref`` is an explicit fixed reference luminosity (Lenstool's ``mag0``), *not* the sample max. The + exponents are fixed at the Faber-Jackson values (L ∝ sigma^4 gives sigma ∝ L^(1/4); constant M/L gives + r_cut ∝ L^(1/2)) — only the normalization ``sigma_ref``, the fiducial velocity dispersion of a + reference-magnitude galaxy, is fitted. Our member luminosities are normalised to the BCG's F160W flux, + so ``L_ref = 1.0`` anchors the relation to the BCG itself. **1 free parameter total for the whole tier — independent of the number of members.** - **Host dark matter halo:** a standalone ``Galaxy`` carrying an ``NFWMCRLudlowSph`` halo with @@ -283,7 +285,7 @@ __Model__ The model is composed below in four blocks: main-tier loop, host halo, source-tier loop, scaling-tier -loop (defining the shared ``b0_ref`` normalization once outside the loop). The four blocks are then +loop (defining the shared ``sigma_ref`` normalization once outside the loop). The four blocks are then bundled into a single ``af.Collection`` model that the analysis will receive. """ redshift_lens = 0.308 @@ -295,11 +297,16 @@ galaxy_models = al.galaxy_af_models_from_csv_tables(mass_table, point_table) -# Main Lens Galaxies: free dPIE ra / rs / b0; centre stays fixed at the CSV value. +# Main Lens Galaxies: free dPIE sigma / r_core / r_cut; centre and redshifts stay +# fixed at the CSV values, and the cosmology constants H0 / Om0 are pinned (they are +# model *constants*, not parameters to sample — if left unset they would inherit the +# config's default priors and float). for name in ("lens_0", "lens_1"): - galaxy_models[name].mass.ra = af.UniformPrior(lower_limit=1.0, upper_limit=15.0) - galaxy_models[name].mass.rs = af.UniformPrior(lower_limit=5.0, upper_limit=40.0) - galaxy_models[name].mass.b0 = af.UniformPrior(lower_limit=0.1, upper_limit=10.0) + galaxy_models[name].mass.sigma = af.UniformPrior(lower_limit=50.0, upper_limit=600.0) + galaxy_models[name].mass.r_core = af.UniformPrior(lower_limit=1.0, upper_limit=15.0) + galaxy_models[name].mass.r_cut = af.UniformPrior(lower_limit=5.0, upper_limit=40.0) + galaxy_models[name].mass.H0 = 67.66 + galaxy_models[name].mass.Om0 = 0.30966 # Host Halo: free mass_at_200; centre + redshift_object + redshift_source fixed. galaxy_models["host_halo"].dark.mass_at_200 = af.LogUniformPrior( @@ -319,19 +326,21 @@ mean=float(np.mean(positions[:, 1])), sigma=3.0 ) -# Scaling Tier (reference-anchored: b0_ref is the single shared free parameter, the -# lens strength of a galaxy at the reference magnitude; per-member ra, rs and b0 -# derive from it with the exponent fixed at the Faber-Jackson value 0.5 — the -# Lenstool convention. The reference luminosity is an EXPLICIT FIXED constant -# (Lenstool's "mag0"); our member luminosities are normalised to the BCG's F160W -# flux, so L_ref = 1.0 anchors the relation to the BCG). +# Scaling Tier (reference-anchored: sigma_ref is the single shared free parameter, +# the fiducial velocity dispersion of a galaxy at the reference magnitude, in km/s; +# per-member sigma, r_core and r_cut derive from it with the exponents fixed at the +# Faber-Jackson values (sigma ∝ L^0.25, radii ∝ L^0.5) — the Lenstool potfile +# convention. The reference luminosity is an EXPLICIT FIXED constant (Lenstool's +# "mag0"); our member luminosities are normalised to the BCG's F160W flux, so +# L_ref = 1.0 anchors the relation to the BCG). -scaling_b0_ref = af.UniformPrior(lower_limit=0.0, upper_limit=1.0) -scaling_exponent = 0.5 +scaling_sigma_ref = af.UniformPrior(lower_limit=0.0, upper_limit=300.0) +scaling_sigma_exponent = 0.25 +scaling_radius_exponent = 0.5 reference_luminosity = 1.0 -scaling_ra_ref_fixed = 0.158 -scaling_rs_ref_fixed = 15.8 +scaling_r_core_ref_fixed = 0.158 +scaling_r_cut_ref_fixed = 15.8 scaling_galaxies_list = [] for centre, luminosity in zip( @@ -341,9 +350,13 @@ mass = af.Model(al.mp.dPIEMassSph) mass.centre = tuple(centre) - mass.ra = scaling_ra_ref_fixed * luminosity_ratio**scaling_exponent - mass.rs = scaling_rs_ref_fixed * luminosity_ratio**scaling_exponent - mass.b0 = scaling_b0_ref * luminosity_ratio**scaling_exponent + mass.sigma = scaling_sigma_ref * luminosity_ratio**scaling_sigma_exponent + mass.r_core = scaling_r_core_ref_fixed * luminosity_ratio**scaling_radius_exponent + mass.r_cut = scaling_r_cut_ref_fixed * luminosity_ratio**scaling_radius_exponent + mass.redshift_object = redshift_lens + mass.redshift_source = max(source_redshifts) + mass.H0 = 67.66 + mass.Om0 = 0.30966 scaling_galaxies_list.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass)) diff --git a/scripts/group/features/linear_light_profiles/slam.py b/scripts/group/features/linear_light_profiles/slam.py index 8b66843aa..b41bf7977 100644 --- a/scripts/group/features/linear_light_profiles/slam.py +++ b/scripts/group/features/linear_light_profiles/slam.py @@ -257,9 +257,21 @@ def source_lp_1( extra_galaxies = af.Collection(extra_mass_models) if extra_mass_models else None - # --- scaling lens galaxy models (light fixed, shared luminosity scaling relation) --- - scaling_factor = af.UniformPrior(lower_limit=0.0, upper_limit=0.5) - scaling_relation = af.UniformPrior(lower_limit=0.0, upper_limit=2.0) + # --- scaling lens galaxy models (light fixed, reference-anchored scaling relation) --- + # Lenstool convention: einstein_radius = einstein_radius_ref * (L / L_ref)^0.5, exponent + # FIXED at the Faber-Jackson value; the single free parameter einstein_radius_ref is the + # Einstein radius of a galaxy as bright as the reference. The reference luminosity is the + # brightest MAIN lens galaxy (the BGG) — only its luminosity anchors the relation, its + # (free) mass is NOT coupled to the tier. See scripts/group/slam.py for the full rationale. + main_luminosity_list = [ + abs(tracer.galaxies[i].bulge.luminosity_within_circle_from(radius=10.0)) + / pixel_scale**2 + for i in range(n_main) + ] + reference_luminosity = max(main_luminosity_list) + + einstein_radius_ref = af.UniformPrior(lower_limit=0.0, upper_limit=2.0) + scaling_exponent = 0.5 scaling_mass_models = [] for i in range(n_scaling): @@ -273,7 +285,10 @@ def source_lp_1( abs(galaxy_with_intensity.bulge.luminosity_within_circle_from(radius=10.0)) / pixel_scale**2 ) - mass.einstein_radius = scaling_factor * total_luminosity**scaling_relation + mass.einstein_radius = ( + einstein_radius_ref + * (total_luminosity / reference_luminosity) ** scaling_exponent + ) scaling_mass_models.append( af.Model( @@ -720,9 +735,18 @@ def mass_total( extra_galaxies = af.Collection(extra_mass_models) if extra_mass_models else None - # --- scaling galaxies: fixed light, free shared scaling relation --- - scaling_factor = af.UniformPrior(lower_limit=0.0, upper_limit=0.5) - scaling_relation = af.UniformPrior(lower_limit=0.0, upper_limit=2.0) + # --- scaling galaxies: fixed light, reference-anchored scaling relation --- + # Same Lenstool convention as source_lp_1: free einstein_radius_ref, exponent fixed + # at 0.5, reference luminosity = the brightest main lens (BGG). + main_luminosity_list = [ + abs(tracer.galaxies[i].bulge.luminosity_within_circle_from(radius=10.0)) + / pixel_scale**2 + for i in range(n_lenses) + ] + reference_luminosity = max(main_luminosity_list) + + einstein_radius_ref = af.UniformPrior(lower_limit=0.0, upper_limit=2.0) + scaling_exponent = 0.5 scaling_mass_models = [] for i in range(n_scaling): @@ -736,7 +760,10 @@ def mass_total( abs(galaxy_with_intensity.bulge.luminosity_within_circle_from(radius=10.0)) / pixel_scale**2 ) - mass.einstein_radius = scaling_factor * total_luminosity**scaling_relation + mass.einstein_radius = ( + einstein_radius_ref + * (total_luminosity / reference_luminosity) ** scaling_exponent + ) scaling_mass_models.append( af.Model( diff --git a/scripts/group/features/scaling_relation/modeling.py b/scripts/group/features/scaling_relation/modeling.py index a8bc5de63..0b4a658c3 100644 --- a/scripts/group/features/scaling_relation/modeling.py +++ b/scripts/group/features/scaling_relation/modeling.py @@ -255,9 +255,10 @@ normalization-slope degeneracy. Only luminosity ratios enter, so the luminosity units are irrelevant; magnitude catalogues convert via ``L / L_ref = 10 ** (0.4 * (m_ref - m))``. -The dPIE-profile cluster-scale analogue — which also scales the core and truncation radii (``ra, rs ∝ L^0.5``, -mirroring Lenstool's r_core / r_cut scaling) — is ``scripts/cluster/modeling.py``. To free the exponent as a -systematics test, replace the fixed value with e.g. ``af.UniformPrior(lower_limit=0.0, upper_limit=1.0)``. +The dPIE-profile cluster-scale analogue — expressed directly in Lenstool's native parameters +(``sigma ∝ L^0.25``, ``r_core, r_cut ∝ L^0.5``, matching Lenstool's potfile scalings; equivalent to +einstein_radius ∝ L^0.5 since the lens strength goes as sigma²) — is ``scripts/cluster/modeling.py``. To free the +exponent as a systematics test, replace the fixed value with e.g. ``af.UniformPrior(lower_limit=0.0, upper_limit=1.0)``. Adding more scaling galaxies (e.g. by lengthening the centres + luminosity lists) does not add any free parameters to the model. diff --git a/scripts/group/slam.py b/scripts/group/slam.py index 5c7559f5b..77e4f1b41 100644 --- a/scripts/group/slam.py +++ b/scripts/group/slam.py @@ -54,12 +54,18 @@ Scaling galaxies are a larger ensemble of companions whose masses are constrained through a shared luminosity-to-mass scaling relation rather than being individually free. They each carry a free MGE -light profile, but their Einstein radii follow: - - einstein_radius = scaling_factor * total_luminosity^scaling_relation - -where `scaling_factor` and `scaling_relation` are two shared free parameters whose priors are -`UniformPrior(0, 0.5)` and `UniformPrior(0, 2)` respectively. This reduces the number of mass +light profile, but their Einstein radii follow the reference-anchored convention used by Lenstool +and essentially every published group- and cluster-scale analysis (Limousin et al. 2005; +Eliasdottir et al. 2007; Bergamini et al. 2019): + + einstein_radius = einstein_radius_ref * (luminosity / reference_luminosity)^0.5 + +The exponent is FIXED at the Faber-Jackson value (einstein_radius ∝ sigma² and sigma ∝ L^(1/4) give +einstein_radius ∝ L^(1/2)), avoiding the normalization-slope degeneracy of fitting both. The tier's +single shared free parameter is `einstein_radius_ref` — the Einstein radius of a galaxy as bright as +the reference — with the reference luminosity anchored to the brightest main lens galaxy (the BGG), +whose luminosity is fixed by the preceding light fit. Only the BGG's *luminosity* anchors the +relation; its (free) mass is deliberately not coupled to the tier. This reduces the number of mass parameters considerably when many companion galaxies are present. Scaling-galaxy models are stored in `model.scaling_galaxies`. @@ -212,9 +218,10 @@ def source_lp_0( Multiple main-lens galaxies each get an `Isothermal` mass; only `lens_0` carries an `ExternalShear`. Extra-galaxy Einstein radii are bounded by a luminosity-derived prior -(`min(5 * 0.5 * L^0.6, 5.0)`). Scaling galaxies share two free parameters, -`scaling_factor` and `scaling_relation`, so their masses follow -`einstein_radius = scaling_factor * luminosity^scaling_relation`. +(`min(5 * 0.5 * L^0.6, 5.0)`). Scaling galaxies share one free parameter, +`einstein_radius_ref`, so their masses follow the reference-anchored Faber-Jackson relation +`einstein_radius = einstein_radius_ref * (luminosity / reference_luminosity)^0.5`, with the +reference luminosity fixed to the brightest main lens (BGG). """ @@ -310,9 +317,30 @@ def source_lp_1( extra_galaxies = af.Collection(extra_mass_models) if extra_mass_models else None - # --- scaling lens galaxy models (light fixed, shared luminosity scaling relation) --- - scaling_factor = af.UniformPrior(lower_limit=0.0, upper_limit=0.5) - scaling_relation = af.UniformPrior(lower_limit=0.0, upper_limit=2.0) + # --- scaling lens galaxy models (light fixed, reference-anchored scaling relation) --- + # Lenstool convention: einstein_radius = einstein_radius_ref * (L / L_ref)^0.5 with the + # exponent FIXED at the Faber-Jackson value (einstein_radius ∝ sigma² and sigma ∝ L^(1/4) + # give einstein_radius ∝ L^(1/2)). The tier's single free parameter is + # einstein_radius_ref — the Einstein radius of a galaxy as bright as the reference — a + # physically interpretable normalization with an easy-to-motivate prior range, avoiding + # the normalization-slope degeneracy of fitting both. The reference luminosity is the + # brightest MAIN lens galaxy (the BGG) — an explicit anchor fixed by the preceding light + # fit and invariant to which galaxies are in the tier. Only the BGG's *luminosity* is the + # anchor; its (free) mass is deliberately NOT coupled to the tier. + main_luminosity_list = [ + np.sum( + [ + 2 * np.pi * g.sigma**2 / g.axis_ratio() * g.intensity + for g in tracer.galaxies[i].bulge.profile_list + ] + ) + / pixel_scale**2 + for i in range(n_main) + ] + reference_luminosity = max(main_luminosity_list) + + einstein_radius_ref = af.UniformPrior(lower_limit=0.0, upper_limit=2.0) + scaling_exponent = 0.5 scaling_mass_models = [] for i in range(n_scaling): @@ -327,7 +355,10 @@ def source_lp_1( for g in tracer.galaxies[n_main + n_extra + i].bulge.profile_list ] total_luminosity = np.sum(luminosity_per_gaussian_list) / pixel_scale**2 - mass.einstein_radius = scaling_factor * total_luminosity**scaling_relation + mass.einstein_radius = ( + einstein_radius_ref + * (total_luminosity / reference_luminosity) ** scaling_exponent + ) scaling_mass_models.append( af.Model( @@ -787,9 +818,25 @@ def mass_total( extra_galaxies = af.Collection(extra_mass_models) if extra_mass_models else None - # --- scaling galaxies: fixed light, free shared scaling relation --- - scaling_factor = af.UniformPrior(lower_limit=0.0, upper_limit=0.5) - scaling_relation = af.UniformPrior(lower_limit=0.0, upper_limit=2.0) + # --- scaling galaxies: fixed light, reference-anchored scaling relation --- + # Same Lenstool convention as source_lp_1: free einstein_radius_ref, exponent fixed + # at 0.5, reference luminosity = the brightest main lens (BGG) from the light[1] fit. + main_luminosity_list = [ + np.sum( + [ + 2 * np.pi * g.sigma**2 / g.axis_ratio() * g.intensity + for g in getattr( + light_result.instance.galaxies, f"lens_{i}" + ).bulge.profile_list + ] + ) + / pixel_scale**2 + for i in range(n_lenses) + ] + reference_luminosity = max(main_luminosity_list) + + einstein_radius_ref = af.UniformPrior(lower_limit=0.0, upper_limit=2.0) + scaling_exponent = 0.5 scaling_mass_models = [] for i in range(n_scaling): @@ -804,7 +851,10 @@ def mass_total( for g in light_scaling.bulge.profile_list ] total_luminosity = np.sum(luminosity_per_gaussian_list) / pixel_scale**2 - mass.einstein_radius = scaling_factor * total_luminosity**scaling_relation + mass.einstein_radius = ( + einstein_radius_ref + * (total_luminosity / reference_luminosity) ** scaling_exponent + ) scaling_mass_models.append( af.Model( diff --git a/scripts/guides/profiles/mass.py b/scripts/guides/profiles/mass.py index 796073bc8..2f481a4da 100644 --- a/scripts/guides/profiles/mass.py +++ b/scripts/guides/profiles/mass.py @@ -123,11 +123,16 @@ power_law_broken = al.mp.PowerLawBroken() power_law_broken_sph = al.mp.PowerLawBrokenSph() -# Pseudo-isothermal family (mass and potential parameterisations) -# Note: dPIEMass with default ell_comps=(0,0) triggers a divide-by-zero in the complex-plane -# formula; we use a small ellipticity here so the survey constructions succeed cleanly. -dpie_mass = al.mp.dPIEMass(ell_comps=(0.05, 0.0)) +# Pseudo-isothermal family. The default dPIEMass / dPIEMassSph are parameterized in +# Lenstool's native convention (ellipticity, angle_pos, sigma = fiducial v_disp in km/s, +# r_core, r_cut, plus the redshifts entering the D_LS/D_S normalization) — the same numbers +# that appear in published cluster papers' results tables. The internal (ra, rs, b0) +# parameterization is the non-standard dPIEMassB0 / dPIEMassB0Sph. +dpie_mass = al.mp.dPIEMass(ellipticity=0.1) dpie_mass_sph = al.mp.dPIEMassSph() +dpie_mass_b0 = al.mp.dPIEMassB0(ell_comps=(0.05, 0.0)) +# Note: PIEMass with ell_comps=(0,0) triggers a divide-by-zero in the complex-plane +# formula; we use a small ellipticity here so the survey constructions succeed cleanly. pie_mass = al.mp.PIEMass(ell_comps=(0.05, 0.0)) dpie_potential = al.mp.dPIEPotential() dpie_potential_sph = al.mp.dPIEPotentialSph() @@ -562,9 +567,11 @@ aplt.plot_array( array=al.mp.dPIEMass( centre=(0.0, 0.0), - ell_comps=al.convert.ell_comps_from(axis_ratio=0.9, angle=45.0), - ra=0.1, - b0=0.5, + ellipticity=0.1, + angle_pos=45.0, + sigma=200.0, + r_core=0.1, + r_cut=20.0, ).convergence_2d_from(grid=grid), title="dPIEMass Convergence", ) diff --git a/scripts/weak/features/strong_lensing/a2744.py b/scripts/weak/features/strong_lensing/a2744.py index 5e94f8761..c17c80e04 100644 --- a/scripts/weak/features/strong_lensing/a2744.py +++ b/scripts/weak/features/strong_lensing/a2744.py @@ -150,9 +150,11 @@ galaxy_models = al.galaxy_af_models_from_csv_tables(mass_table, point_table) for name in ("lens_0", "lens_1"): - galaxy_models[name].mass.ra = af.UniformPrior(lower_limit=1.0, upper_limit=15.0) - galaxy_models[name].mass.rs = af.UniformPrior(lower_limit=5.0, upper_limit=40.0) - galaxy_models[name].mass.b0 = af.UniformPrior(lower_limit=0.1, upper_limit=10.0) + galaxy_models[name].mass.sigma = af.UniformPrior(lower_limit=50.0, upper_limit=600.0) + galaxy_models[name].mass.r_core = af.UniformPrior(lower_limit=1.0, upper_limit=15.0) + galaxy_models[name].mass.r_cut = af.UniformPrior(lower_limit=5.0, upper_limit=40.0) + galaxy_models[name].mass.H0 = 67.66 + galaxy_models[name].mass.Om0 = 0.30966 galaxy_models["host_halo"].dark.mass_at_200 = af.LogUniformPrior( lower_limit=10**14.5, upper_limit=10**16.0 @@ -168,11 +170,13 @@ mean=float(np.mean(positions[:, 1])), sigma=3.0 ) -scaling_b0_ref = af.UniformPrior(lower_limit=0.0, upper_limit=1.0) -scaling_exponent = 0.5 +scaling_sigma_ref = af.UniformPrior(lower_limit=0.0, upper_limit=300.0) +scaling_sigma_exponent = 0.25 +scaling_radius_exponent = 0.5 reference_luminosity = 1.0 -scaling_ra_ref_fixed = 0.158 -scaling_rs_ref_fixed = 15.8 +scaling_r_core_ref_fixed = 0.158 +scaling_r_cut_ref_fixed = 15.8 +source_redshift_max = max(float(d.redshift) for d in dataset_list) scaling_galaxies_list = [] for centre, luminosity in zip( @@ -182,9 +186,13 @@ mass = af.Model(al.mp.dPIEMassSph) mass.centre = tuple(centre) - mass.ra = scaling_ra_ref_fixed * luminosity_ratio**scaling_exponent - mass.rs = scaling_rs_ref_fixed * luminosity_ratio**scaling_exponent - mass.b0 = scaling_b0_ref * luminosity_ratio**scaling_exponent + mass.sigma = scaling_sigma_ref * luminosity_ratio**scaling_sigma_exponent + mass.r_core = scaling_r_core_ref_fixed * luminosity_ratio**scaling_radius_exponent + mass.r_cut = scaling_r_cut_ref_fixed * luminosity_ratio**scaling_radius_exponent + mass.redshift_object = redshift_lens + mass.redshift_source = source_redshift_max + mass.H0 = 67.66 + mass.Om0 = 0.30966 scaling_galaxies_list.append(af.Model(al.Galaxy, redshift=redshift_lens, mass=mass)) diff --git a/workspace_index.json b/workspace_index.json index 9f75acef9..f84482698 100644 --- a/workspace_index.json +++ b/workspace_index.json @@ -202,6 +202,7 @@ "autolens_workspace/start_here.py", "cluster/modeling.py", "csv_api.py", + "prep.py", "scripts/cluster/csv_api.py", "scripts/cluster/modeling.py", "scripts/cluster/simulator.py", @@ -1471,6 +1472,21 @@ "summary": "Over sampling is a numerical technique where the images of light profiles and galaxies are evaluated on a higher resolution grid than the image data to ensure the calculation is accurate.", "title": "Chaining: Over Sample" }, + { + "contents": [ + "Simulate", + "Smooth Model Fit", + "Joint Fit", + "Dkappa Map", + "Iterative Fit", + "Evidence Sampling" + ], + "cross_refs": [], + "notebook": "notebooks/guides/advanced/potential_correction.ipynb", + "path": "scripts/guides/advanced/potential_correction.py", + "summary": "Parametric lens mass models (e.g. the `PowerLaw` used throughout the workspace) are smooth: their convergence and deflection angles vary gradually across the image plane. Real lens galaxies are not perfectly smooth. They contain substructure \u2014 dark matter subhaloes, globular clusters, satellite galaxies \u2014 and their large-scale mass distributions can depart from ellipsoidal symmetry in ways no simple parametric profile captures.", + "title": "Potential Corrections (Gravitational Imaging)" + }, { "contents": [ "Conventions", @@ -3895,6 +3911,40 @@ "summary": "Source science focuses on studying the highly magnified properties of the background lensed source galaxy (or galaxies).", "title": "Pixelization: Source Science" }, + { + "contents": [ + "Simplifications", + "Prerequisites", + "Dataset", + "Arc Mask", + "Smooth Starting Model", + "Dpsi Mesh", + "Dpsi Gradient Operator", + "Source Gradients", + "PSF Blur Matrix", + "Dpsi Mapping Matrix", + "Source Inversion Blocks", + "Joint System", + "Solve", + "Dkappa Map", + "Evidence Terms", + "Verification", + "JAX / xp", + "Wrap Up" + ], + "cross_refs": [ + "autolens_workspace_test/scripts/jax_likelihood_functions/imaging/potential_correction.py", + "features/pixelization/likelihood_function.py", + "guides/advanced/potential_correction.ipynb", + "guides/advanced/potential_correction.py", + "imaging/features/pixelization/likelihood_function.ipynb", + "imaging/likelihood_function.ipynb" + ], + "notebook": "notebooks/imaging/features/potential_correction/likelihood_function.ipynb", + "path": "scripts/imaging/features/potential_correction/likelihood_function.py", + "summary": "This script provides a step-by-step guide of the **PyAutoLens** potential-correction `log_likelihood_function` (the Bayesian evidence of `al.pc.FitDpsiSrcImaging`), which jointly inverts `Imaging` data for a pixelized source and pixelized corrections $\\delta\\psi$ to the lensing potential.", + "title": "__Log Likelihood Function: Potential Correction (Gravitational Imaging)__" + }, { "contents": [ "Prerequisites", @@ -5017,6 +5067,49 @@ "summary": "Source science focuses on studying the highly magnified properties of the background lensed source galaxy (or galaxies).", "title": "Pixelization: Source Reconstruction" }, + { + "contents": [ + "Simplifications", + "Prerequisites", + "Dataset", + "Source Inversion Blocks", + "Dpsi Mesh", + "Correction Response", + "Joint Response", + "Dense Route", + "Sparse (w-tilde) Route", + "Solve", + "Dkappa Map", + "Evidence Terms", + "Verification", + "JAX / xp", + "Wrap Up" + ], + "cross_refs": [ + "autolens_workspace_test/scripts/jax_likelihood_functions/interferometer/potential_correction.py", + "imaging/features/potential_correction/likelihood_function.ipynb", + "interferometer/features/pixelization/likelihood_function.ipynb", + "interferometer/features/pixelization/likelihood_function.py", + "interferometer/features/potential_correction/start_here.ipynb", + "start_here.py" + ], + "notebook": "notebooks/interferometer/features/potential_correction/likelihood_function.ipynb", + "path": "scripts/interferometer/features/potential_correction/likelihood_function.py", + "summary": "This script provides a step-by-step guide of the **PyAutoLens** visibility-space potential-correction `log_likelihood_function` (the Bayesian evidence of `al.pc.FitDpsiSrcInterferometer`), which jointly inverts `Interferometer` data for a pixelized source and pixelized corrections $\\delta\\psi$ to the lensing potential.", + "title": "__Log Likelihood Function: Potential Correction (Interferometer)__" + }, + { + "contents": [], + "cross_refs": [ + "guides/advanced/potential_correction.ipynb", + "interferometer/start_here.ipynb", + "interferometer/start_here.py" + ], + "notebook": "notebooks/interferometer/features/potential_correction/start_here.ipynb", + "path": "scripts/interferometer/features/potential_correction/start_here.py", + "summary": "This example performs gravitational imaging on interferometer data: pixelized corrections $\\delta\\psi$ to the lensing potential are reconstructed jointly with the pixelized source, revealing mass structure (e.g. dark subhaloes) the smooth lens model omits, as a convergence-correction map $\\delta\\kappa = \\frac{1}{2}\\nabla^2\\delta\\psi$.", + "title": "Feature: Potential Correction (Gravitational Imaging)" + }, { "contents": [ "SLaM Pipelines", @@ -5299,6 +5392,24 @@ "summary": "This script simulates `Interferometer` data of a 'galaxy-scale' strong lens where:", "title": "Simulator: Mutli Interferometer" }, + { + "contents": [ + "Dataset", + "Quasar Image Masking", + "Model", + "Analysis Factor & Factor Graph", + "Search & Model-Fit", + "Result" + ], + "cross_refs": [ + "point_source/features/deblending.py", + "point_source/start_here.py" + ], + "notebook": "notebooks/multi/features/imaging_and_point_source/modeling.ipynb", + "path": "scripts/multi/features/imaging_and_point_source/modeling.py", + "summary": "This script jointly fits CCD imaging of a lensed quasar's **extended arcs** and its **point-source observables** (image positions and time delays) with a **single lens mass model**, using PyAutoFit's factor-graph API to sample the joint likelihood with one non-linear search.", + "title": "Modeling: Combined Imaging + Point Source" + }, { "contents": [ "Model", @@ -5824,7 +5935,27 @@ "notebook": "notebooks/point_source/start_here.ipynb", "path": "scripts/point_source/start_here.py", "summary": "Strong gravitational lenses often have point sources (e.g. quasars) that are being lensed, appearing as two or four distinct point-like images. These lenses are particularly useful for measuring cosmological parameters like the Hubble constant, and for studying the small-scale properties of dark matter.", - "title": "Start Here: Imaging" + "title": "Start Here: Point Source" + }, + { + "contents": [ + "Dataset (Strong)", + "Dataset (Weak)", + "Model", + "Analysis Factors & Factor Graph", + "Search & Model-Fit", + "Result" + ], + "cross_refs": [ + "cluster/start_here.py", + "modeling.py", + "weak/features/strong_lensing/modeling.py", + "weak/start_here.py" + ], + "notebook": "notebooks/weak/features/strong_lensing/a2744.ipynb", + "path": "scripts/weak/features/strong_lensing/a2744.py", + "summary": "This script constrains a **single cluster-scale mass model** of **Abell 2744** (\"Pandora's Cluster\", z = 0.308) using both its **strong-lensing** multiple images and its **weak-lensing** shear catalogue \u2014 real data on both sides, fitted jointly with one non-linear search via PyAutoFit's factor-graph API.", + "title": "Real Data: Combined Strong + Weak Lensing of Abell 2744" }, { "contents": [ @@ -5978,5 +6109,34 @@ "path": "scripts/weak/simulator.py", "summary": "This script simulates a weak gravitational lensing shear catalogue. Unlike the imaging simulator (which produces a 2D image of the lensed source) the weak-lensing simulator produces a *catalogue* of (gamma_2, gamma_1) shear measurements at the (y, x) positions of a population of background source galaxies.", "title": "Simulator: Weak Lensing" + }, + { + "contents": [ + "JAX", + "Google Colab Setup", + "Imports", + "Catalogue Download", + "Catalogue Load & Projection", + "Quality Cuts", + "Weak Dataset", + "Mass Map", + "Model", + "Model Fit", + "Result", + "Model Your Own Catalogue", + "Simulator", + "Wrap Up" + ], + "cross_refs": [ + "/weak/likelihood_function.py", + "/weak/modeling.py", + "/weak/real_data/a2744.py", + "autolens_workspace/start_here.py", + "scripts/weak/simulator.py" + ], + "notebook": "notebooks/weak/start_here.ipynb", + "path": "scripts/weak/start_here.py", + "summary": "Gravitational lensing does not stop at the spectacular arcs and multiple images of strong lensing. Every galaxy behind a massive structure is slightly sheared by its gravity \u2014 a percent-level distortion invisible in any single galaxy but measurable statistically across thousands. This is **weak lensing**, and it traces mass on scales far beyond the strong-lensing core.", + "title": "Start Here: Weak Lensing" } ]