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

def test__convergence__correct_values(self):
        # eta = 1.0
        # kappa = 0.5 * 1.0 ** 1.0
        isothermal = mp.SphericalIsothermal(centre=(0.0, 0.0), einstein_radius=2.0)
        assert isothermal.convergence_from_grid(grid=np.array([[0.0, 1.0]])) == pytest.approx(0.5 * 2.0, 1e-3)

        isothermal = mp.EllipticalIsothermal(centre=(0.0, 0.0), axis_ratio=1.0, phi=0.0, einstein_radius=1.0)
        assert isothermal.convergence_from_grid(grid=np.array([[0.0, 1.0]])) == pytest.approx(0.5, 1e-3)

        isothermal = mp.EllipticalIsothermal(centre=(0.0, 0.0), axis_ratio=1.0, phi=0.0, einstein_radius=2.0)
        assert isothermal.convergence_from_grid(grid=np.array([[0.0, 1.0]])) == pytest.approx(0.5 * 2.0, 1e-3)

        isothermal = mp.EllipticalIsothermal(centre=(0.0, 0.0), axis_ratio=0.5, phi=0.0, einstein_radius=1.0)
        assert isothermal.convergence_from_grid(grid=np.array([[0.0, 1.0]])) == pytest.approx(0.66666, 1e-3)

def test__mass_within_ellipse__compare_to_grid__uses_conversion_factor(self):

        sie = mp.EllipticalIsothermal(einstein_radius=2.0, axis_ratio=0.5, phi=0.0)

        integral_radius = 0.5
        dimensionless_mass_tot = 0.0

        xs = np.linspace(-1.0, 1.0, 40)
        ys = np.linspace(-1.0, 1.0, 40)

        edge = xs[1] - xs[0]
        area = edge ** 2

        for x in xs:
            for y in ys:

                eta = sie.grid_to_elliptical_radii(np.array([[x, y]]))

                if eta < integral_radius:

def test__compare_radial_critical_curves_from_magnification_and_lamda_t__reg_grid_two_component_galaxy(
            self
        ):
            grid = grids.Grid.from_shape_pixel_scale_and_sub_grid_size(
                shape=(100, 100), pixel_scale=0.05, sub_grid_size=1
            )

            g0 = g.Galaxy(
                redshift=0.5,
                mass_profile=mp.EllipticalIsothermal(
                    centre=(0.0, 0.0), einstein_radius=1.4, axis_ratio=0.7, phi=40.0
                ),
            )

            g1 = g.Galaxy(
                redshift=0.5,
                mass_profile=mp.SphericalIsothermal(
                    centre=(1.0, 1.0), einstein_radius=2.0
                ),
            )

            plane = pl.Plane(galaxies=[g0, g1], redshift=None)

            critical_curve_radial_from_magnification = critical_curve_via_magnification_from_plane_and_grid(
                plane=plane, grid=grid
            )[

def test__compare_radial_caustic_from_magnification_and_lambda_t__two_component_galaxy(
            self
        ):
            grid = grids.Grid.from_shape_pixel_scale_and_sub_grid_size(
                shape=(60, 60), pixel_scale=0.5, sub_grid_size=2
            )

            g0 = g.Galaxy(
                redshift=0.5,
                mass_profile=mp.EllipticalIsothermal(
                    centre=(0.0, 0.0), einstein_radius=1.4, axis_ratio=0.7, phi=40.0
                ),
            )

            g1 = g.Galaxy(
                redshift=0.5,
                mass_profile=mp.SphericalIsothermal(
                    centre=(1.0, 1.0), einstein_radius=2.0
                ),
            )

            plane = pl.Plane(galaxies=[g0, g1], redshift=None)

            caustic_radial_from_magnification = caustics_via_magnification_from_plane_and_grid(
                plane=plane, grid=grid
            )[

print('Grid')
print(lens_data.grid_stack.regular)

# The image, noise-map and grids are masked using the mask and mapped to 1D arrays for fast calcuations.
print(lens_data.image.shape) # This is the original 2D image
print(lens_data.image_1d.shape)
print(lens_data.noise_map_1d.shape)
print(lens_data.grid_stack.regular.shape)

# To fit an image, we need to create an image-plane image using a tracer.
# Lets use the same tracer we simulated the ccd data with (thus, our fit should be 'perfect').

# Its worth noting that below, we use the lens_data's grid-stack to setup the tracer. This ensures that our image-plane
# image will be the same resolution and alignment as our image-data, as well as being masked appropriately.

lens_galaxy = g.Galaxy(mass=mp.EllipticalIsothermal(centre=(0.0, 0.0), einstein_radius=1.6, axis_ratio=0.7, phi=45.0))
source_galaxy = g.Galaxy(light=lp.EllipticalSersic(centre=(0.1, 0.1), axis_ratio=0.8, phi=45.0,
                                                        intensity=1.0, effective_radius=1.0, sersic_index=2.5))
tracer = ray_tracing.TracerImageSourcePlanes(lens_galaxies=[lens_galaxy], source_galaxies=[source_galaxy],
                                             image_plane_grid_stack=lens_data.grid_stack)
ray_tracing_plotters.plot_image_plane_image(tracer=tracer)

# To fit the image, we pass the lens data and tracer to the fitting module. This performs the following:

# 1) Blurs the tracer's image-plane image with the lens data's PSF, ensuring that the telescope optics are
#    accounted for by the fit. This creates the fit's 'model_image'.

# 2) Computes the difference between this model_image and the observed image-data, creating the fit's 'residual_map'.

# 3) Divides the residuals by the noise-map and squaring each value, creating the fit's 'chi_squared_map'.

# 4) Sums up these chi-squared values and converts them to a 'likelihood', which quantities how good the tracer's fit

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.05, pixel_scale=0.05)
    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.8, phi=45.0,
                                                        einstein_radius=0.8))

    source_galaxy = g.Galaxy(light=lp.EllipticalSersic(centre=(0.0, 0.0), axis_ratio=0.5, phi=90.0, intensity=0.03,
                                                       effective_radius=0.3, sersic_index=3.0))
    tracer = ray_tracing.TracerImageSourcePlanes(lens_galaxies=[lens_galaxy], source_galaxies=[source_galaxy],
                                                 image_plane_grid_stack=image_plane_grid_stack)

    ccd_simulated = 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)

    return ccd_simulated

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.05, pixel_scale=0.05)

    image_plane_grid_stack = grids.GridStack.grid_stack_for_simulation(shape=(180, 180), pixel_scale=0.05,
                                                                       psf_shape=(11, 11))

    lens_galaxy = g.Galaxy(mass=mp.EllipticalIsothermal(centre=(0.0, 0.0), axis_ratio=0.8, phi=135.0,
                                                        einstein_radius=1.6))

    source_galaxy = g.Galaxy(light=lp.EllipticalSersic(centre=(0.1, 0.1), axis_ratio=0.8, phi=90.0, intensity=0.2,
                                                       effective_radius=0.3, sersic_index=1.0))

    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.05,
                                exposure_time=300.0, psf=psf, background_sky_level=0.1, add_noise=True)

def pass_priors(self, previous_results):

            phase_1_results = previous_results[0]
            phase_2_results = previous_results[1]

            # We're going to link the centres of the light profiles computed above to the centre of the lens galaxy
            # mass-profiles in this phase. Because the centres of the mass profiles were fixed in phases 1 and 2,
            # linking them using the 'variable' attribute means that they stay constant (which for now, is what we want).

            self.lens_galaxies.left_lens.mass.centre_0 = phase_1_results.variable.left_lens.light.centre_0
            self.lens_galaxies.left_lens.mass.centre_1 = phase_1_results.variable.left_lens.light.centre_1

            self.lens_galaxies.right_lens.mass.centre_0 = phase_2_results.variable.right_lens.light.centre_0
            self.lens_galaxies.right_lens.mass.centre_1 = phase_2_results.variable.right_lens.light.centre_1

    phase3 = LensSubtractedPhase(lens_galaxies=dict(left_lens=gm.GalaxyModel(mass=mp.EllipticalIsothermal),
                                                    right_lens=gm.GalaxyModel(mass=mp.EllipticalIsothermal)),
                                 source_galaxies=dict(source=gm.GalaxyModel(light=lp.EllipticalExponential)),
                                 optimizer_class=nl.MultiNest,
                                 phase_name=pipeline_path + '/phase_3_fit_sources')

    phase3.optimizer.const_efficiency_mode = True
    phase3.optimizer.n_live_points = 50
    phase3.optimizer.sampling_efficiency = 0.5

    ### PHASE 4 ###

    # In phase 4, we'll fit both lens galaxy's light and mass profiles, as well as the source-galaxy, simultaneously.

    class FitAllPhase(ph.LensSourcePlanePhase):

        def pass_priors(self, previous_results):

phase2.optimizer.const_efficiency_mode = True
    phase2.optimizer.n_live_points = 40
    phase2.optimizer.sampling_efficiency = 0.5

    # Now lets do the same again, but with 3 source galaxy components.

    class X3SourcePhase(ph.LensSourcePlanePhase):

        def pass_priors(self, previous_results):

            self.lens_galaxies.lens = previous_results[1].variable.lens
            self.source_galaxies.source.light_0 = previous_results[1].variable.source.light_0
            self.source_galaxies.source.light_1 = previous_results[1].variable.source.light_1

    phase3 = X3SourcePhase(lens_galaxies=dict(lens=gm.GalaxyModel(mass=mp.EllipticalIsothermal)),
                           source_galaxies=dict(source=gm.GalaxyModel(light_0=lp.EllipticalExponential,
                                                                      light_1=lp.EllipticalSersic,
                                                                      light_2=lp.EllipticalSersic)),
                           optimizer_class=nl.MultiNest, phase_name=pipeline_path + '/phase_3_x3_source')

    phase3.optimizer.const_efficiency_mode = True
    phase3.optimizer.n_live_points = 50
    phase3.optimizer.sampling_efficiency = 0.5

    # And one more for luck!

    class X4SourcePhase(ph.LensSourcePlanePhase):

        def pass_priors(self, previous_results):

            self.lens_galaxies.lens = previous_results[2].variable.lens

# 5) A tracer_without_subhalo can make an regular-plane + source-plane strong lens system.
# 6) The Universe's cosmology can be input into this tracer_without_subhalo to convert units to physical values.
# 7) The tracer_without_subhalo's regular-plane regular can be used to simulate strong lens ccd observed on a real telescope.
# 8) That this regular can be fitted, so to as quantify how well a model strong lens system represents the observed regular.

# In this summary, we'll consider how flexible the tools PyAutoLens gives you are to study every aspect of a strong
# lens system. Lets get a 'fit_normal' to a strong lens, by setting up an regular, masks, tracer_without_subhalo, etc.

path = 'path/to/AutoLens/howtolens/chapter_1_introduction' # Unfortunately, in a Jupyter notebook you have to manually specify the path to PyAutoLens and this tutorial.
path = '/home/jammy/PyCharm/Projects/AutoLens/workspace/howtolens/chapter_1_introduction'
image = im.load_ccd_data_from_fits(image_path=path + '/datas/regular.fits',
                                   noise_map_path=path+'/datas/noise_maps.fits',
                                   psf_path=path + '/datas/psf.fits', pixel_scale=0.1)
mask = ma.Mask.circular(shape=image.shape, pixel_scale=image.pixel_scale, radius_arcsec=3.0)
lensing_image = li.LensData(ccd_data=image, mask=mask)
lens_galaxy = g.Galaxy(mass=mp.EllipticalIsothermal(centre=(0.0, 0.0), einstein_radius=1.6, axis_ratio=0.7, phi=45.0))
source_galaxy = g.Galaxy(bulge=lp.EllipticalSersic(centre=(0.1, 0.1), axis_ratio=0.8, phi=45.0,
                                                  intensity=1.0, effective_radius=1.0, sersic_index=4.0),
                         disk=lp.EllipticalSersic(centre=(0.1, 0.1), axis_ratio=0.8, phi=45.0,
                                                  intensity=1.0, effective_radius=1.0, sersic_index=1.0))
tracer = ray_tracing.TracerImageSourcePlanes(lens_galaxies=[lens_galaxy], source_galaxies=[source_galaxy],
                                             image_plane_grid_stack=[lensing_image.grid_stack])
fit = lens_fit.fit_lens_image_with_tracer(lens_image=lensing_image, tracer=tracer)

# The fit_normal contains our tracer_without_subhalo, which contains our planes, which contain our grid_stacks and galaxies, which contain our
# profiles:
print(fit)
print()
print(fit.tracer)
print()
print(fit.tracer.image_plane)
print()