How to use the autofit.optimize.non_linear.MultiNest function in autofit

def test__phase_can_receive_list_of_galaxy_models(self):
        phase = ph.LensPlanePhase(lens_galaxies=[gm.GalaxyModel(sersic=lp.EllipticalSersic,
                                                                sis=mp.SphericalIsothermal,
                                                                variable_redshift=True),
                                                 gm.GalaxyModel(sis=mp.SphericalIsothermal,
                                                                variable_redshift=True)],
                                  optimizer_class=non_linear.MultiNest, phase_name='test_phase')

        instance = phase.optimizer.variable.instance_from_physical_vector(
            [0.2, 0.21, 0.22, 0.23, 0.24, 0.25, 0.8, 0.1, 0.2, 0.3,
             0.4, 0.9, 0.5, 0.7, 0.8])

        assert instance.lens_galaxies[0].sersic.centre[0] == 0.2
        assert instance.lens_galaxies[0].sis.centre[0] == 0.1
        assert instance.lens_galaxies[0].sis.centre[1] == 0.2
        assert instance.lens_galaxies[0].sis.einstein_radius == 0.3
        assert instance.lens_galaxies[0].redshift == 0.4
        assert instance.lens_galaxies[1].sis.centre[0] == 0.9
        assert instance.lens_galaxies[1].sis.centre[1] == 0.5
        assert instance.lens_galaxies[1].sis.einstein_radius == 0.7
        assert instance.lens_galaxies[1].redshift == 0.8

        class LensPlanePhase2(ph.LensPlanePhase):

# We can also customize the source galaxy - lets say we believe it is compact and limit its effective radius

        self.source_galaxies.source_galaxy.light.effective_radius = \
            model_mapper.UniformPrior(lower_limit=0.0, upper_limit=0.3)

# (If the use of a python 'class' here, or some of the python code doesn't appear clear, don't worry about it. This
# code is using a number of Python's object-oriented features. In general, I expect that'll you'll simply copy the code
# above and use it as a template).

# We can now create this custom phase like we did a normal phase before. When we run the phase, the pass_prior function
# will be called automatically and thus change the priors as we specified above. If you look at the 'model.info'
# file in the output of the non-linear search, you'll see that the priors have indeed been changed.
custom_phase = CustomPhase(lens_galaxies=dict(lens_galaxy=lens_galaxy_model),
                           source_galaxies=dict(source_galaxy=source_galaxy_model),
                           optimizer_class=non_linear.MultiNest,
                           phase_name='2_custom_priors')

print('MultiNest has begun running - checkout the workspace/howtolens/chapter_2_lens_modeling/output/2_custom_priors'
      'folder for live output of the results, images and lens model.'
      'This Jupyter notebook cell with progress once MultiNest has completed - this could take some time!')
results_custom = custom_phase.run(data=ccd_data)
lens_fit_plotters.plot_fit_subplot(fit=results_custom.most_likely_fit)
print('MultiNest has finished run - you may now continue the notebook.')

# In phase 1, we will fit the lens galaxy's mass and one source galaxy, where we:

    # 1) Set our priors on the lens galaxy (y,x) centre such that we assume the image is centred around the lens galaxy.

    class LensSourceX1Phase(ph.LensSourcePlanePhase):

        def pass_priors(self, previous_results):

            self.lens_galaxies.lens.mass.centre_0 = prior.GaussianPrior(mean=0.0, sigma=0.1)
            self.lens_galaxies.lens.mass.centre_1 = prior.GaussianPrior(mean=0.0, sigma=0.1)

    phase1 = LensSourceX1Phase(lens_galaxies=dict(lens=gm.GalaxyModel(mass=mp.EllipticalIsothermal,
                                                                      shear=mp.ExternalShear)),
                               source_galaxies=dict(source_0=gm.GalaxyModel(light=lp.EllipticalSersic)),
                               mask_function=mask_function, optimizer_class=nl.MultiNest,
                               phase_name=pipeline_path + '/phase_1_x1_source')

    # You'll see these lines throughout all of the example pipelines. They are used to make MultiNest sample the \
    # non-linear parameter space faster (if you haven't already, checkout 'tutorial_7_multinest_black_magic' in
    # 'howtolens/chapter_2_lens_modeling'.

    # Fitting the lens galaxy and source galaxy from uninitialized priors often risks MultiNest getting stuck in a
    # local maxima, especially for the image in this example which actually has two source galaxies. Therefore, whilst
    # I will continue to use constant efficiency mode to ensure fast run time, I've upped the number of live points
    # and decreased the sampling efficiency from the usual values to ensure the non-linear search is robust.

    phase1.optimizer.const_efficiency_mode = True
    phase1.optimizer.n_live_points = 80
    phase1.optimizer.sampling_efficiency = 0.2

    ### PHASE 2 ###

return ccd_simulated

# Simulate the image and set it up.
ccd_data = simulate()
ccd_plotters.plot_ccd_subplot(ccd_data=ccd_data)

# Lets first run the phase without black magic, which is performed using the standard lines of code we're used to.

# A word of warning, this phase takes >1 hour to run... so if you get bored, skip the run cell below continue to the
# phase with black magic (afterall, the whole point of this tutorial is how slow MultiNest can run, so theres no harm
# if the slow run speed bores you to tears :P).

phase_normal = ph.LensSourcePlanePhase(lens_galaxies=dict(lens=gm.GalaxyModel(light=lp.EllipticalSersic,
                                                                                   mass=mp.EllipticalIsothermal)),
                                            source_galaxies=dict(source=gm.GalaxyModel(light=lp.EllipticalSersic)),
                                            optimizer_class=nl.MultiNest, phase_name='7_no_black_magic')

# We're going to use the time module to time how long each MultiNest run takes. However, if you resume the MultiNest
# run from a previous job, this time won't be accurate. Fortunately, if you look in the folder
# 'howtolens/chapter_2_lens_modeling/output/7_no_black_magic') you'll find a file 'run_time' which gives the overall
# run-time including any resumes.
start = time.time()

# Lets run the phase - the run-time will be output to the output/7_multinest_black_magic/

print('MultiNest has begun running - checkout the workspace/howtolens/chapter_2_lens_modeling/output/7_multinest_black_magic'
      ' folder for live output of the results, images and lens model.'
      ' This Jupyter notebook cell with progress once MultiNest has completed - this could take some time!')
phase_normal_results = phase_normal.run(data=ccd_data)
print('MultiNest has finished run - you may now continue the notebook.')

# Lets check we get a reasonably good model and fit to the data.

# When we plot it, the lens light's is clealy visible in the centre of the image
ccd_plotters.plot_ccd_subplot(ccd_data=ccd_data)

# Now lets fit it using a phase, noting that indeed the galaxy-model corresponds to the one above.

# Because we now have 18 non-linear parameters, the non-linear search takes a lot longer to run. On my laptop, this
# phase took around an hour, which is a bit too long for you to wait if you want to go through these tutorials quickly.
# Therefore, as discussed before, I've included the results of this non-linear search already, allowing you to
# go through the tutorial as if you had actually run them.

# Nevertheless, you could try running it yourself (maybe over your lunch break?). All you need to do is change the
# phase_name below, maybe to something like 'howtolens/3_realism_and_complexity_rerun'
phase = ph.LensSourcePlanePhase(lens_galaxies=dict(lens_galaxy=gm.GalaxyModel(light=lp.EllipticalSersic,
                                                                                mass=mp.EllipticalIsothermal)),
                                source_galaxies=dict(source_galaxy=gm.GalaxyModel(light=lp.EllipticalExponential)),
                                optimizer_class=nl.MultiNest,
                                phase_name='3_realism_and_complexity')

# Lets run the phase.
print('MultiNest has begun running - checkout the workspace/howtolens/chapter_2_lens_modeling/output/3_realism_and_complexity'
      'folder for live output of the results, images and lens model.'
      'This Jupyter notebook cell with progress once MultiNest has completed - this could take some time!')
results = phase.run(data=ccd_data)
print('MultiNest has finished run - you may now continue the notebook.')

# And lets look at the image.
lens_fit_plotters.plot_fit_subplot(fit=results.most_likely_fit, should_plot_mask=True, zoom_around_mask=True)

# Uh-oh. That image didn't look very good, did it? If we compare our inferred parameters (look at the
# 'workspace/howtolens/chapter_2_lens_modeling/output/3_realism_and_complexity' folder to the actual
# values (in the simulate function) you'll see that we have, indeed, fitted the wrong model.

# In the pass priors function, we can 'pair' any two parameters by setting them equal to one another. This
        # removes the parameter on the left-hand side of the pairing from the lens model, such that is always assumes
        # the same value as the parameter on the right-hand side.
        self.lens_galaxies.lens.mass.centre_0 = self.lens_galaxies.lens.light.centre_0

        # Now, the mass-profile's x coordinate will only use the x coordinate of the light profile. Lets do this with
        # the remaining geometric parameters of the light and mass profiles
        self.lens_galaxies.lens.mass.centre_1 = self.lens_galaxies.lens.light.centre_1
        self.lens_galaxies.lens.mass.axis_ratio = self.lens_galaxies.lens.light.axis_ratio
        self.lens_galaxies.lens.mass.phi = self.lens_galaxies.lens.light.phi

# Again, we create this phase and run it. The non-linear search has a less complex parameter space to seach, and thus
light_traces_mass_phase = LightTracesMassPhase(lens_galaxies=dict(lens=gm.GalaxyModel(light=lp.EllipticalSersic,
                                                                                       mass=mp.EllipticalIsothermal)),
                                      source_galaxies=dict(source=gm.GalaxyModel(light=lp.EllipticalExponential)),
                                      optimizer_class=nl.MultiNest,
                                               phase_name='4_light_traces_mass')
print('MultiNest has begun running - checkout the workspace/howtolens/chapter_2_lens_modeling/output/4_dealing_with_failure'
      'folder for live output of the results, images and lens model.'
      'This Jupyter notebook cell with progress once MultiNest has completed - this could take some time!')
light_traces_mass_phase_result = light_traces_mass_phase.run(data=ccd_data)
print('MultiNest has finished run - you may now continue the notebook.')
lens_fit_plotters.plot_fit_subplot(fit=light_traces_mass_phase_result.most_likely_fit)

# The results look pretty good. Our source galaxy fits the data pretty well, and we've clearly inferred a model that
# looks similar to the one above. However, inspection of the residuals shows that the fit wasn't quite as good as the
# custom-phase above.

# It turns out that when I simulated this image, light didn't perfectly trace mass. The light-profile's axis-ratio was
# 0.9, whereas the mass-profiles was 0.8. The quality of the fit has suffered as a result, and the likelihood we've
# inferred is lower.
#

self.lens_galaxies.los_0.mass.centre_1 = previous_results[0].constant.los_0.light.centre_1
            self.lens_galaxies.los_1.mass.centre_0 = previous_results[0].constant.los_1.light.centre_0
            self.lens_galaxies.los_1.mass.centre_1 = previous_results[0].constant.los_1.light.centre_1
            self.lens_galaxies.los_2.mass.centre_0 = previous_results[0].constant.los_2.light.centre_0
            self.lens_galaxies.los_2.mass.centre_1 = previous_results[0].constant.los_2.light.centre_1

            self.source_galaxies.source = previous_results[2].variable.source

    phase4 = SingleLensPlanePhase(lens_galaxies=dict(lens=gm.GalaxyModel(mass=mp.EllipticalIsothermal),
                                                     los_0=gm.GalaxyModel(mass=mp.SphericalIsothermal),
                                                     los_1=gm.GalaxyModel(mass=mp.SphericalIsothermal),
                                                     los_2=gm.GalaxyModel(mass=mp.SphericalIsothermal)),
                                  source_galaxies=dict(source=gm.GalaxyModel(pixelization=pix.AdaptiveMagnification,
                                                                             regularization=reg.Constant)),
                                  use_positions=True,
                                  optimizer_class=nl.MultiNest, phase_name=pipeline_path + '/phase_4_single_plane')

    # Customize MultiNest so it runs fast
    phase4.optimizer.n_live_points = 60
    phase4.optimizer.sampling_efficiency = 0.2
    phase4.optimizer.const_efficiency_mode = True


    # In phase 5, we will fit the foreground subtracted image using a multi-plane ray tracer. This means that the
    # redshift of each line-of-sight galaxy is included as a free parameter (we assume the lens and source redshifts
    # are known).

    class MultiPlanePhase(ph.MultiPlanePhase):

        def modify_image(self, image, previous_results):
            return image - previous_results[0].unmasked_lens_plane_model_image

# 1) Set our priors on the lens galaxy (y,x) centre such that we assume the image is centred around the lens galaxy.
    # 2) Use a circular mask which therefore also includes the source's light.

    # Whereas in the howtolens tutorial I used an anti-annular mask in this phase to remove the source galaxy, here I
    # use the default 3.0"  circular mask. In general, I haven't found the choice of mask to make a big difference,
    # albeit this does depend on how much off the lens galaxy's light the lensed source galaxy's light obstructs.

    class LensPhase(ph.LensPlanePhase):

        def pass_priors(self, previous_results):

            self.lens_galaxies.lens.light.centre_0 = mm.GaussianPrior(mean=0.0, sigma=0.1)
            self.lens_galaxies.lens.light.centre_1 = mm.GaussianPrior(mean=0.0, sigma=0.1)

    phase1 = LensPhase(lens_galaxies=dict(lens=gm.GalaxyModel(light=lp.EllipticalSersic)),
                       optimizer_class=nl.MultiNest, phase_name=pipeline_path + '/phase_1_lens_light_only')

    # You'll see these lines throughout all of the example pipelines. They are used to make MultiNest sample the \
    # non-linear parameter space faster (if you haven't already, checkout 'tutorial_7_multinest_black_magic' in
    # 'howtolens/chapter_2_lens_modeling'.

    phase1.optimizer.const_efficiency_mode = True
    phase1.optimizer.n_live_points = 30
    phase1.optimizer.sampling_efficiency = 0.3

    ### PHASE 2 ###

    # In phase 2, we will fit the lens galaxy's mass and source galaxy's light, where we:

    # 1) Use a lens-subtracted image generated by subtracting model lens galaxy image from phase 1.
    # 2) Use a circular annular mask to only include the source-galaxy.
    # 3) Initialize the priors on the centre of the lens galaxy's mass-profile by linking them to those inferred for \

# included! They're the pixels the source is actually located. Therefore, we're going to use an 'anti-annular
    # mask', where the inner and outer radii are the regions we omit from the analysis. This means we need to specify
    # a third ring of the masks, even further out, such that data at these exterior edges of the image are also masked.

    # We can change a mask using the 'mask_function', which basically returns the new masks we want to use (you can
    # actually use on phases by themselves like in the previous chapter).

    def mask_function(image):
        return msk.Mask.circular_anti_annular(shape=image.shape, pixel_scale=image.pixel_scale, inner_radius_arcsec=0.5,
                                              outer_radius_arcsec=1.6, outer_radius_2_arcsec=2.5)

    # We next create the phase, using the same notation we learnt before (but noting the masks function is passed to
    # this phase ensuring the anti-annular masks above is used).

    phase1 = phase.LensPlanePhase(lens_galaxies=dict(lens=gm.GalaxyModel(light=lp.EllipticalSersic)),
                                  optimizer_class=nl.MultiNest, mask_function=mask_function,
                                  phase_name=pipeline_path + '/phase_1_lens_light_only')

    # We'll use the MultiNest black magic we covered in tutorial 7 of chapter 2 to get this phase to run fast.

    phase1.optimizer.const_efficiency_mode = True
    phase1.optimizer.n_live_points = 30
    phase1.optimizer.sampling_efficiency = 0.5

    ### PHASE 2 ###

    # In phase 2, we fit the source galaxy's light. Thus, we want to make 2 changes from the previous phase

    # 1) We want to fit the lens subtracted image calculated in phase 1, instead of the observed image.
    # 2) We want to mask the central regions of this image, where there are residuals due to the lens light subtraction.

    # We can use the mask function again, to modify the mask to an annulus. We'll use the same ring radii as before.

class X4SourcePhase(ph.LensSourcePlanePhase):

        def pass_priors(self, previous_results):

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

    phase4 = X4SourcePhase(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,
                                                                      light_3=lp.EllipticalSersic)),
                           optimizer_class=nl.MultiNest, phase_name=pipeline_path + '/phase_4_x4_source')

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

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