How to use the autolens.model.profiles.mass_profiles function in autolens

def test__constructor(self):
        nfw = mp.EllipticalNFW(centre=(0.7, 1.0), axis_ratio=0.7, phi=60.0, kappa_s=2.0,
                               scale_radius=10.0)

        assert nfw.centre == (0.7, 1.0)
        assert nfw.axis_ratio == 0.7
        assert nfw.phi == 60.0
        assert nfw.kappa_s == 2.0
        assert nfw.scale_radius == 10.0

        nfw = mp.SphericalNFW(centre=(0.7, 1.0), kappa_s=2.0, scale_radius=10.0)

        assert nfw.centre == (0.7, 1.0)
        assert nfw.axis_ratio == 1.0
        assert nfw.phi == 0.0
        assert nfw.kappa_s == 2.0
        assert nfw.scale_radius == 10.0

def test__compare_to_cored_power_law(self):
        power_law = mp.EllipticalCoredIsothermal(centre=(0.0, 0.0), axis_ratio=0.5, phi=45.0, einstein_radius=1.0,
                                                 core_radius=0.1)
        cored_power_law = mp.EllipticalCoredPowerLaw(centre=(0.0, 0.0), axis_ratio=0.5, phi=45.0,
                                                     einstein_radius=1.0, slope=2.0, core_radius=0.1)

        assert power_law.potential_from_grid(grid) == pytest.approx(cored_power_law.potential_from_grid(grid), 1e-3)
        assert power_law.potential_from_grid(grid) == pytest.approx(cored_power_law.potential_from_grid(grid), 1e-3)
        assert power_law.deflections_from_grid(grid) == pytest.approx(cored_power_law.deflections_from_grid(grid), 1e-3)
        assert power_law.deflections_from_grid(grid) == pytest.approx(cored_power_law.deflections_from_grid(grid), 1e-3)

mass_to_light_ratio=1.0)
        assert dev.convergence_from_grid(grid=np.array([[1.0, 0.0]])) == pytest.approx(5.6697, 1e-3)

        dev = mp.EllipticalDevVaucouleurs(axis_ratio=0.5, phi=90.0, intensity=2.0, effective_radius=3.0,
                                          mass_to_light_ratio=1.0)
        assert dev.convergence_from_grid(grid=np.array([[0.0, 1.0]])) == pytest.approx(7.4455, 1e-3)

        dev = mp.EllipticalDevVaucouleurs(axis_ratio=0.5, phi=90.0, intensity=4.0, effective_radius=3.0,
                                          mass_to_light_ratio=1.0)
        assert dev.convergence_from_grid(grid=np.array([[0.0, 1.0]])) == pytest.approx(2.0 * 7.4455, 1e-3)

        dev = mp.EllipticalDevVaucouleurs(axis_ratio=0.5, phi=90.0, intensity=2.0, effective_radius=3.0,
                                          mass_to_light_ratio=2.0)
        assert dev.convergence_from_grid(grid=np.array([[0.0, 1.0]])) == pytest.approx(2.0 * 7.4455, 1e-3)

        sersic = mp.SphericalDevVaucouleurs(centre=(0.0, 0.0), intensity=1.0, effective_radius=0.6,
                                            mass_to_light_ratio=1.0)
        assert sersic.convergence_from_grid(grid=np.array([[0.0, 1.0]])) == pytest.approx(0.351797, 1e-3)

def test__potential_correct_values(self):
        gnfw = mp.SphericalGeneralizedNFW(centre=(0.0, 0.0), kappa_s=1.0, inner_slope=0.5, scale_radius=8.0)
        assert gnfw.potential_from_grid(grid=np.array([[0.1625, 0.1875]])) == pytest.approx(0.00920, 1e-3)

        gnfw = mp.SphericalGeneralizedNFW(centre=(0.0, 0.0), kappa_s=1.0, inner_slope=1.5, scale_radius=8.0)
        assert gnfw.potential_from_grid(grid=np.array([[0.1625, 0.1875]])) == pytest.approx(0.17448, 1e-3)

        gnfw = mp.EllipticalGeneralizedNFW(centre=(1.0, 1.0), kappa_s=5.0, axis_ratio=0.5,
                                           phi=100.0, inner_slope=1.0, scale_radius=10.0)
        assert gnfw.potential_from_grid(grid=np.array([[2.0, 2.0]])) == pytest.approx(2.4718, 1e-4)

def test__deflections__correct_values(self):

        power_law = mp.SphericalPowerLaw(centre=(0.2, 0.2), einstein_radius=1.0, slope=2.0)
        defls = power_law.deflections_from_grid(grid=np.array([[0.1875, 0.1625]]))
        assert defls[0, 0] == pytest.approx(-0.31622, 1e-3)
        assert defls[0, 1] == pytest.approx(-0.94868, 1e-3)

        power_law = mp.SphericalPowerLaw(centre=(0.2, 0.2), einstein_radius=1.0, slope=2.5)
        defls = power_law.deflections_from_grid(grid=np.array([[0.1875, 0.1625]]))
        assert defls[0, 0] == pytest.approx(-1.59054, 1e-3)
        assert defls[0, 1] == pytest.approx(-4.77162, 1e-3)

        power_law = mp.SphericalPowerLaw(centre=(0.2, 0.2), einstein_radius=1.0, slope=1.5)
        defls = power_law.deflections_from_grid(grid=np.array([[0.1875, 0.1625]]))
        assert defls[0, 0] == pytest.approx(-0.06287, 1e-3)
        assert defls[0, 1] == pytest.approx(-0.18861, 1e-3)

        power_law = mp.EllipticalPowerLaw(centre=(0, 0), axis_ratio=0.5, phi=0.0, einstein_radius=1.0, slope=2.5)
        defls = power_law.deflections_from_grid(grid=np.array([[0.1625, 0.1625]]))
        assert defls[0, 0] == pytest.approx(1.29641, 1e-3)
        assert defls[0, 1] == pytest.approx(0.99629, 1e-3)

        power_law = mp.EllipticalPowerLaw(centre=(0, 0), axis_ratio=0.5, phi=0.0, einstein_radius=1.0, slope=1.5)
        defls = power_law.deflections_from_grid(grid=np.array([[0.1625, 0.1625]]))
        assert defls[0, 0] == pytest.approx(0.48036, 1e-3)
        assert defls[0, 1] == pytest.approx(0.26729, 1e-3)

        power_law = mp.EllipticalPowerLaw(centre=(-0.7, 0.5), axis_ratio=0.7, phi=60.0, einstein_radius=1.3,
                                          slope=1.9)

def test__deflections__change_geometry(self):
        
        point_mass_0 = mp.PointMass(centre=(0.0, 0.0))
        point_mass_1 = mp.PointMass(centre=(1.0, 1.0))
        defls_0 = point_mass_0.deflections_from_grid(grid=np.array([[1.0, 1.0]]))
        defls_1 = point_mass_1.deflections_from_grid(grid=np.array([[0.0, 0.0]]))
        assert defls_0[0, 0] == pytest.approx(-defls_1[0, 0], 1e-5)
        assert defls_0[0, 1] == pytest.approx(-defls_1[0, 1], 1e-5)

        point_mass_0 = mp.PointMass(centre=(0.0, 0.0))
        point_mass_1 = mp.PointMass(centre=(0.0, 0.0))
        defls_0 = point_mass_0.deflections_from_grid(grid=np.array([[1.0, 0.0]]))
        defls_1 = point_mass_1.deflections_from_grid(grid=np.array([[0.0, 1.0]]))
        assert defls_0[0, 0] == pytest.approx(defls_1[0, 1], 1e-5)
        assert defls_0[0, 1] == pytest.approx(defls_1[0, 0], 1e-5)

def test__x1_lens_galaxy__other_childs_of_power_law_all_work(self, grid_stack):

            g0 = g.Galaxy(mass=mp.EllipticalIsothermal(einstein_radius=1.0), redshift=1.0)

            tracer = ray_tracing.TracerImageSourcePlanes(lens_galaxies=[g0],
                                                         source_galaxies=[g.Galaxy(redshift=2.0)],
                                                         image_plane_grid_stack=grid_stack,
                                                         cosmology=cosmo.Planck15)

            g0_mass = g0.mass_within_circle(radius=1.0, conversion_factor=tracer.critical_density_arcsec)

            assert tracer.einstein_mass_of_lens_galaxy == g0_mass

            g0 = g.Galaxy(mass=mp.SphericalCoredIsothermal(einstein_radius=1.0, core_radius=0.1), redshift=1.0)

            tracer = ray_tracing.TracerImageSourcePlanes(lens_galaxies=[g0],
                                                         source_galaxies=[g.Galaxy(redshift=2.0)],
                                                         image_plane_grid_stack=grid_stack,
                                                         cosmology=cosmo.Planck15)

# In phase 3, we will fit simultaneously the lens and source galaxies, where we:

    # 1) Initialize the lens's light, mass, shear and source's light using the results of phases 1 and 2.

    class LensSourcePhase(ph.LensSourcePlanePhase):

        def pass_priors(self, previous_results):

            self.lens_galaxies.lens.light = previous_results[0].variable.lens.light
            self.lens_galaxies.lens.mass = previous_results[1].variable.lens.mass
            self.lens_galaxies.lens.shear = previous_results[1].variable.lens.shear
            self.source_galaxies.source = previous_results[1].variable.source

    phase3 = LensSourcePhase(lens_galaxies=dict(lens=gm.GalaxyModel(light=lp.EllipticalSersic,
                                                                    mass=mp.EllipticalIsothermal,
                                                                    shear=mp.ExternalShear)),
                             source_galaxies=dict(source=gm.GalaxyModel(light=lp.EllipticalSersic)),
                             optimizer_class=nl.MultiNest, phase_name=pipeline_path + '/phase_3_both')

    phase3.optimizer.const_efficiency_mode = True
    phase3.optimizer.n_live_points = 75
    phase3.optimizer.sampling_efficiency = 0.3

    return pipeline.PipelineImaging(pipeline_path, phase1, phase2, phase3)

# - If two or more source galaxies are at the same redshift in the source-plane, their light is ray-traced in the
#   same way. Therefore, when determining their lensed images, we can sum the lensed images of each galaxy's light-profiles.

# So, lets do it - lets use the 'plane' module in AutoLens to create a strong lensing system like the one pictured
# above. For simplicity, we'll assume 1 lens galaxy and 1 source galaxy.

# As always, we need grids, where our grids are the coordinates we'll 'trace' from the image-plane to the source-plane
# in the lensing configuration above. Our grid-stack is therefore no longer just a 'grid-stack', but the grid-stack
# representing our image-plane coordinates. Thus, lets name as such.
image_plane_grid_stack = grids.GridStack.from_shape_pixel_scale_and_sub_grid_size(shape=(100, 100), pixel_scale=0.05,
                                                                                  sub_grid_size=2)

# Whereas before we called our galaxy's things like 'galaxy_with_light_profile', lets now refer to them by their role
# in lensing, e.g. 'lens_galaxy' and 'source_galaxy'.
mass_profile = mass_profiles.SphericalIsothermal(centre=(0.0,  0.0), einstein_radius=1.6)
lens_galaxy = galaxy.Galaxy(mass=mass_profile)
light_profile = light_profiles.SphericalSersic(centre=(0.0, 0.0), intensity=1.0, effective_radius=1.0, sersic_index=1.0)
source_galaxy = galaxy.Galaxy(light=light_profile)

# Lets setup our image-plane. This plane takes the lens galaxy we made above and the grid-stack of
# image-plane coordinates.
image_plane = plane.Plane(galaxies=[lens_galaxy], grid_stack=image_plane_grid_stack)

# Up to now, we've kept our galaxies and grids separate, and passed the grid to a galaxy object to compute
# its quantities (e.g. to compute light-profile intensities, we'd write galaxy.intensities_from_grid(grid=grid)).
#
# Plane's combine the galaxies and grids into one object, thus once we've setup a plane there is no longer any need
# have to pass it a grid to compute its quantities. Furthermore, once these quantities are in a plane, they are
# automatically mapped back to their original 2D grid-arrays.
print('deflection-angles of planes regular-grid pixel 1:')
print(image_plane.deflections_y[0,0])

def simulate():

    from autolens.data.array import grids
    from autolens.model.galaxy import galaxy as g
    from autolens.lens import ray_tracing

    psf = ccd.PSF.simulate_as_gaussian(shape=(11, 11), sigma=0.1, pixel_scale=0.1)

    image_plane_grid_stack = grids.GridStack.grid_stack_for_simulation(shape=(130, 130), pixel_scale=0.1,
                                                                       psf_shape=(11, 11))

    lens_galaxy = g.Galaxy(light=lp.EllipticalSersic(centre=(0.0, 0.0), axis_ratio=0.9, phi=45.0, intensity=0.04,
                                                     effective_radius=0.5, sersic_index=3.5),
                           mass=mp.EllipticalIsothermal(centre=(0.0, 0.0), axis_ratio=0.9, phi=45.0,
                                                        einstein_radius=1.0),
                           shear=mp.ExternalShear(magnitude=0.05, phi=90.0))

    source_galaxy = g.Galaxy(light=lp.SphericalExponential(centre=(0.0, 0.0), intensity=0.2, effective_radius=0.2))
    tracer = ray_tracing.TracerImageSourcePlanes(lens_galaxies=[lens_galaxy], source_galaxies=[source_galaxy],
                                                 image_plane_grid_stack=image_plane_grid_stack)

    return ccd.CCDData.simulate(array=tracer.image_plane_image_for_simulation, pixel_scale=0.1,
                                exposure_time=300.0, psf=psf, background_sky_level=0.1, add_noise=True)