diff --git a/config/non_linear/mle.yaml b/config/non_linear/mle.yaml index 365b911fb..2cc331c7f 100644 --- a/config/non_linear/mle.yaml +++ b/config/non_linear/mle.yaml @@ -2,53 +2,9 @@ # **PyAutoFit** supports the following maximum likelihood estimator (MLE) algorithms: -# - PySwarms: https://github.com/ljvmiranda921/pyswarms / https://pyswarms.readthedocs.io/en/latest/index.html - # Settings in the [search], [run] and [options] entries are specific to each nested algorithm and should be # determined by consulting that method's own readthedocs. -PySwarmsGlobal: - run: - iters: 2000 - search: - cognitive: 0.5 - ftol: -.inf - inertia: 0.9 - n_particles: 50 - social: 0.3 - initialize: # The method used to generate where walkers are initialized in parameter space {prior | ball}. - method: ball # priors: samples are initialized by randomly drawing from each parameter's prior. ball: samples are initialized by randomly drawing unit values from a narrow uniform distribution. - ball_lower_limit: 0.49 # The lower limit of the uniform distribution unit values are drawn from when initializing walkers using the ball method. - ball_upper_limit: 0.51 # The upper limit of the uniform distribution unit values are drawn from when initializing walkers using the ball method. - parallel: - number_of_cores: 1 # The number of cores the search is parallelized over by default, using Python multiprocessing. - printing: - silence: false # If True, the default print output of the non-linear search is silcened and not printed by the Python interpreter. - iterations_per_full_update: 500 # Non-linear search iterations between every full update, which outputs all visuals and result fits (e.g. model.result, search.summary), this exits the search and can be slow. - iterations_per_quick_update: 500 # Non-linear search iterations between every quick update, which just displays the maximum likelihood model fit. - remove_state_files_at_end: true # Whether to remove the savestate of the seach (e.g. the Emcee hdf5 file) at the end to save hard-disk space (results are still stored as PyAutoFit pickles and loadable). -PySwarmsLocal: - run: - iters: 2000 - search: - cognitive: 0.5 - ftol: -.inf - inertia: 0.9 - minkowski_p_norm: 2 - n_particles: 50 - number_of_k_neighbors: 3 - social: 0.3 - initialize: # The method used to generate where walkers are initialized in parameter space {prior | ball}. - method: ball # priors: samples are initialized by randomly drawing from each parameter's prior. ball: samples are initialized by randomly drawing unit values from a narrow uniform distribution. - ball_lower_limit: 0.49 # The lower limit of the uniform distribution unit values are drawn from when initializing walkers using the ball method. - ball_upper_limit: 0.51 # The upper limit of the uniform distribution unit values are drawn from when initializing walkers using the ball method. - parallel: - number_of_cores: 1 # The number of cores the search is parallelized over by default, using Python multiprocessing. - printing: - silence: false # If True, the default print output of the non-linear search is silcened and not printed by the Python interpreter. - iterations_per_full_update: 500 # Non-linear search iterations between every full update, which outputs all visuals and result fits (e.g. model.result, search.summary), this exits the search and can be slow. - iterations_per_quick_update: 500 # Non-linear search iterations between every quick update, which just displays the maximum likelihood model fit. - remove_state_files_at_end: true # Whether to remove the savestate of the seach (e.g. the Emcee hdf5 file) at the end to save hard-disk space (results are still stored as PyAutoFit pickles and loadable). LBFGS: search: tol: null diff --git a/config/non_linear/nest.yaml b/config/non_linear/nest.yaml index 67a71bd7f..141c3f565 100644 --- a/config/non_linear/nest.yaml +++ b/config/non_linear/nest.yaml @@ -2,7 +2,6 @@ # - Nautilus https://https://github.com/johannesulf/nautilus / https://nautilus-sampler.readthedocs.io/en/stable/index.html # - Dynesty: https://github.com/joshspeagle/dynesty / https://dynesty.readthedocs.io/en/latest/index.html -# - UltraNest: https://github.com/JohannesBuchner/UltraNest / https://johannesbuchner.github.io/UltraNest/readme.html # Settings in the [search] and [run] entries are specific to each nested algorithm and should be determined by # consulting that MCMC method's own readthedocs. Nautilus: @@ -90,51 +89,3 @@ DynestyDynamic: printing: silence: false # If True, the default print output of the non-linear search is silenced and not printed by the Python interpreter. -UltraNest: - search: - draw_multiple: true - ndraw_max: 65536 - ndraw_min: 128 - num_bootstraps: 30 - num_test_samples: 2 - resume: true - run_num: null - storage_backend: hdf5 - vectorized: false - warmstart_max_tau: -1.0 - run: - cluster_num_live_points: 40 - dkl: 0.5 - dlogz: 0.5 - frac_remain: 0.01 - insertion_test_window: 10 - insertion_test_zscore_threshold: 2 - lepsilon: 0.001 - log_interval: null - max_iters: null - max_ncalls: null - max_num_improvement_loops: -1.0 - min_ess: 400 - min_num_live_points: 400 - show_status: true - update_interval_ncall: null - update_interval_volume_fraction: 0.8 - viz_callback: auto - stepsampler: - adaptive_nsteps: false - log: false - max_nsteps: 1000 - nsteps: 25 - region_filter: false - scale: 1.0 - stepsampler_cls: null - initialize: # The method used to generate where walkers are initialized in parameter space {prior}. - method: prior # priors: samples are initialized by randomly drawing from each parameter's prior. - parallel: - number_of_cores: 1 # The number of cores the search is parallelized over by default, using Python multiprocessing. - printing: - silence: false # If True, the default print output of the non-linear search is silenced and not printed by the Python interpreter. - - iterations_per_full_update: 1000 # Non-linear search iterations between every full update, which outputs all visuals and result fits (e.g. model.result, search.summary), this exits the search and can be slow. - iterations_per_quick_update: 500 # Non-linear search iterations between every quick update, which just displays the maximum likelihood model fit. - remove_state_files_at_end: true # Whether to remove the savestate of the seach (e.g. the Emcee hdf5 file) at the end to save hard-disk space (results are still stored as PyAutoFit pickles and loadable). diff --git a/config/visualize/plots_search.yaml b/config/visualize/plots_search.yaml index 0915cfb82..7e61f9536 100644 --- a/config/visualize/plots_search.yaml +++ b/config/visualize/plots_search.yaml @@ -10,15 +10,6 @@ emcee: trajectories: true # Output Emcee trajectories figure during a non-linear search fit? nautilus: cornerplot: true # Output Nautilus cornerplot figure during a non-linear search fit? -pyswarms: - contour: true # Output PySwarms contour figure during a non-linear search fit? - cost_history: true # Output PySwarms cost_history figure during a non-linear search fit? - time_series: true # Output PySwarms time_series figure during a non-linear search fit? - trajectories: true # Output PySwarms trajectories figure during a non-linear search fit? -ultranest: - cornerplot: true # Output Ultranest cornerplot figure during a non-linear search fit? - runplot: true # Output Ultranest runplot figure during a non-linear search fit? - traceplot: true # Output Ultranest traceplot figure during a non-linear search fit? zeus: corner: true # Output Zeus corner figure during a non-linear search fit? likelihood_series: true # Output Zeus likelihood series figure during a non-linear search fit? diff --git a/scripts/guides/modeling/searches.py b/scripts/guides/modeling/searches.py index 2ee35beab..f5e3e9fa8 100644 --- a/scripts/guides/modeling/searches.py +++ b/scripts/guides/modeling/searches.py @@ -26,7 +26,7 @@ **Dynesty:** Dynesty (https://github.com/joshspeagle/dynesty) is a nested sampling algorithm. **Emcee:** Emcee (https://github.com/dfm/emcee) is an ensemble MCMC sampler that is commonly used in Astronomy. **Zeus:** Zeus (https://zeus-mcmc.readthedocs.io/en/latest/) is an ensemble MCMC slice sampler. -**PySwarms:** PySwarms is a particle swarm optimization (PSO) algorithm. +**LBFGS:** LBFGS is a quasi-Newton optimization algorithm from scipy. **Start Point:** For maximum likelihood estimator (MLE) and Markov Chain Monte Carlo (MCMC) non-linear searches. **Search Cookbook:** There are a number of other searches supported by **PyAutoFit** and therefore which can be used. @@ -43,7 +43,6 @@ # %cd $workspace_path # print(f"Working Directory has been set to `{workspace_path}`") -import numpy as np from pathlib import Path import autofit as af import autolens as al @@ -174,34 +173,21 @@ ) """ -__PySwarms__ +__LBFGS__ -PySwarms is a particle swarm optimization (PSO) algorithm. +LBFGS is a quasi-Newton optimization algorithm from scipy. -Information about PySwarms can be found at the following links: +An optimizer only seeks to find the maximum likelihood lens model, unlike MCMC or nested sampling algorithms +like Zeus and Nautilus, which aim to map out parameter space and infer errors on the parameters. Therefore, in +principle, an optimizer like LBFGS should fit a lens model very fast. - - https://github.com/ljvmiranda921/pyswarms - - https://pyswarms.readthedocs.io/en/latest/index.html - - https://pyswarms.readthedocs.io/en/latest/api/pyswarms.single.html#module-pyswarms.single.global_best - -An PSO algorithm only seeks to only find the maximum likelihood lens model, unlike MCMC or nested sampling algorithms -like Zzeus and Nautilus, which aims to map-out parameter space and infer errors on the parameters.Therefore, in -principle, a PSO like PySwarm should fit a lens model very fast. - -In our experience, the parameter spaces fitted by lens models are too complex for `PySwarms` to be used without a lot -of user attention and care. +In our experience, the parameter spaces fitted by lens models are often too complex for optimizers to be used without +careful initialization. """ -search = af.PySwarmsGlobal( +search = af.LBFGS( path_prefix=Path("imaging", "searches"), - name="PySwarmsGlobal", + name="LBFGS", unique_tag="example", - n_particles=30, - iters=300, - cognitive=0.5, - social=0.3, - inertia=0.9, - ftol=-np.inf, - iterations_per_quick_update=1000, ) """ @@ -295,7 +281,7 @@ __Search Cookbook__ There are a number of other searches supported by **PyAutoFit** and therefore which can be used, which are not -explictly documented here. These include Ultranest and LBFGS. +explictly documented here. These include LBFGS. The **PyAutoFit** search cookbook documents all searches that are available, including those not documented here, and provides the code you can easily copy and paste to use these methods. diff --git a/scripts/guides/plot/examples/searches.py b/scripts/guides/plot/examples/searches.py index 64eb46e98..438d4ab65 100644 --- a/scripts/guides/plot/examples/searches.py +++ b/scripts/guides/plot/examples/searches.py @@ -18,7 +18,6 @@ **EmceePlotter:** We now use the MCMC plotting functions to visualize emcee's results. **Search Specific Visualization:** The internal sampler can be used to plot the results of the non-linear search. **ZeusPlotter:** We now use the MCMC plotting functions to visualize Zeus's results. -**PySwarmsPlotter:** We now use the MLE plotting functions to visualize pyswarms's results. **GetDist:** This example illustrates how to plot visualization summarizing the results of model-fit using any. **Parameter Names:** Note that in order to customize the figure, we will use the `samples.model.parameter_names` list. **GetDist Plotter:** To make plots we use a GetDist plotter object, which can be customized to change the appearance of. @@ -629,63 +628,6 @@ except AttributeError: pass # MCMC-specific methods not available for non-MCMC searches -""" -__PySwarmsPlotter__ - -We now use the MLE plotting functions to visualize pyswarms's results. - -The pyswarms readthedocs describes fully all of the methods used below - - - https://pyswarms.readthedocs.io/en/latest/api/pyswarms.utils.plotters.html - -In all the examples below, we use the `kwargs` of this function to pass in any of the input parameters that are -described in the API docs. -""" - -""" -__Search Specific Visualization__ - -The internal sampler can be used to plot the results of the non-linear search. - -We do this using the `search_internal` attribute which contains the sampler in its native form. - -The first time you run a search, the `search_internal` attribute will be available because it is passed ot the -result via memory. - -If you rerun the fit on a completed result, it will not be available in memory, and therefore -will be loaded from the `files/search_internal` folder. The `search_internal` entry of the `output.yaml` must be true -for this to be possible. -""" -search_internal = result.search_internal - -""" -__Plots__ - -The `contour` method shows a 2D projection of the particle trajectories. -""" -try: - from pyswarms.utils import plotters - - plotters.plot_contour( - pos_history=search_internal.pos_history, - canvas=None, - title="Trajectories", - mark=None, - designer=None, - mesher=None, - animator=None, - ) - plt.show() - - plotters.plot_cost_history( - cost_history=search_internal.cost_history, - ax=None, - title="Cost History", - designer=None, - ) - plt.show() -except (AttributeError, ImportError, TypeError): - pass # pyswarms or search_internal unavailable """ __GetDist__ diff --git a/scripts/howtolens/chapter_optional/tutorial_searches.py b/scripts/howtolens/chapter_optional/tutorial_searches.py index db3b148db..25ce892c7 100644 --- a/scripts/howtolens/chapter_optional/tutorial_searches.py +++ b/scripts/howtolens/chapter_optional/tutorial_searches.py @@ -14,7 +14,7 @@ __Contents__ **Nested Sampling:** Lets first perform the model-fit using Nautilus, but look at different parameters that control how. -**Optimizers:** Nested sampling algorithms like Nautilus provides the errors on all of the model parameters, by. +**Optimizers:** There are a class of non-linear searches called optimizers, which seek to optimize the log likelihood. **MCMC:** For users familiar with Markov Chain Monte Carlo (MCMC) non-linear samplers, PyAutoFit supports the. """ @@ -209,94 +209,24 @@ error estimates (perhaps this is our final lens model fit before we publish the results in a paper), these extra iterations are acceptable. -However, we often don't care about the errors. For example, in the previous tutorial when chaining searches, the only +However, we often don't care about the errors. For example, in the previous tutorial when chaining searches, the only result we used from the fit performed in the first search was the maximum log likelihood model, omitting the errors entirely! Its seems wasteful to use a nested sampling algorithm like Nautilus to map out the entirity of parameter -space when we don't use this information! +space when we don't use this information! -There are a class of non-linear searches called `optimizers`, which seek to optimize just one thing, the log -likelihood. They want to find the model that maximizes the log likelihood, with no regard for the errors, thus not -wasting time mapping out in intricate detail every facet of parameter space. Lets see how much faster we can find a -good fit to the lens data using an optimizer. +There are a class of non-linear searches called `optimizers`, which seek to optimize just one thing, the log +likelihood. They want to find the model that maximizes the log likelihood, with no regard for the errors, thus not +wasting time mapping out in intricate detail every facet of parameter space. -we'll use the `Particle Swarm Optimizer` PySwarms. Conceptually this works quite similar to Nautilus, it has a set of -points in parameter space (called `particles`) and it uses their likelihoods to determine where it thinks the higher -likelihood regions of parameter space are. - -Unlike Nautilus, this algorithm requires us to specify how many iterations it should perform to find the global -maxima solutions. Here, an iteration is the number of samples performed by every particle, so the total number of -iterations is n_particles * iters. Lets try a total of 50000 iterations, a factor 10 less than our Nautilus runs above. - -In our experience, pyswarms is ineffective at initializing a lens model and therefore needs a the initial swarm of -particles to surround the highest likelihood lens models. We set this starting point up below by manually inputting -`GaussianPriors` on every parameter, where the centre of these priors is near the true values of the simulated lens data. - -Given this need for a robust starting point, PySwarms is only suited to model-fits where we have this information. It may -therefore be useful when performing lens modeling search chaining (see HowToLens chapter 3). However, even in such -circumstances, we have found that is often unrealible and often infers a local maxima. -""" -lens_bulge = af.Model(al.lp.Sersic) -lens_bulge.centre.centre_0 = af.GaussianPrior(mean=0.0, sigma=0.3) -lens_bulge.centre.centre_1 = af.GaussianPrior(mean=0.0, sigma=0.3) -lens_bulge.ell_comps.ell_comps_0 = af.GaussianPrior(mean=0.0, sigma=0.3) -lens_bulge.ell_comps.ell_comps_1 = af.GaussianPrior(mean=0.0, sigma=0.3) -lens_bulge.intensity = af.GaussianPrior(mean=1.0, sigma=0.3) -lens_bulge.effective_radius = af.GaussianPrior(mean=0.8, sigma=0.2) -lens_bulge.sersic_index = af.GaussianPrior(mean=4.0, sigma=1.0) - -mass = af.Model(al.mp.Isothermal) -mass.centre.centre_0 = af.GaussianPrior(mean=0.0, sigma=0.1) -mass.centre.centre_1 = af.GaussianPrior(mean=0.0, sigma=0.1) -mass.ell_comps.ell_comps_0 = af.GaussianPrior(mean=0.0, sigma=0.3) -mass.ell_comps.ell_comps_1 = af.GaussianPrior(mean=0.0, sigma=0.3) -mass.einstein_radius = af.GaussianPrior(mean=1.4, sigma=0.4) - -shear = af.Model(al.mp.ExternalShear) -shear.gamma_1 = af.GaussianPrior(mean=0.0, sigma=0.1) -shear.gamma_2 = af.GaussianPrior(mean=0.0, sigma=0.1) - -bulge = af.Model(al.lp.Sersic) -bulge.centre.centre_0 = af.GaussianPrior(mean=0.0, sigma=0.3) -bulge.centre.centre_1 = af.GaussianPrior(mean=0.0, sigma=0.3) -bulge.ell_comps.ell_comps_0 = af.GaussianPrior(mean=0.0, sigma=0.3) -bulge.ell_comps.ell_comps_1 = af.GaussianPrior(mean=0.0, sigma=0.3) -bulge.intensity = af.GaussianPrior(mean=0.3, sigma=0.3) -bulge.effective_radius = af.GaussianPrior(mean=0.2, sigma=0.2) -bulge.sersic_index = af.GaussianPrior(mean=1.0, sigma=1.0) - -lens = af.Model(al.Galaxy, redshift=0.5, mass=mass, shear=shear) -source = af.Model(al.Galaxy, redshift=1.0, bulge=bulge) - -model = af.Collection(galaxies=af.Collection(lens=lens, source=source)) - -search = af.PySwarmsLocal( - path_prefix=Path("howtolens", "chapter_optional"), - name="tutorial_searches_pso", - unique_tag=dataset_name, - n_particles=50, - iters=1000, -) - -print( - "The non-linear search has begun running - checkout the workspace/output" - " folder for live output of the results, images and lens model." - " This Jupyter notebook cell with progress once search has completed - this could take some time!" -) - -result_pso = search.fit(model=model, analysis=analysis) - -print("PySwarms has finished run - you may now continue the notebook.") - -aplt.subplot_fit_imaging(fit=result_pso.max_log_likelihood_fit) - -""" -In our experience, the parameter spaces fitted by lens models are too complex for `PySwarms` to be used without a lot -of user attention and care and careful setting up of the initialization priors, as shown above. +PyAutoFit supports the LBFGS optimizer (from scipy), which can be used as an alternative to nested sampling when +only the maximum likelihood model is needed. However, optimizers generally need a good starting point to work well, +and in our experience the parameter spaces fitted by lens models are often too complex for optimizers without careful +setup of initialization priors. __MCMC__ For users familiar with Markov Chain Monte Carlo (MCMC) non-linear samplers, PyAutoFit supports the non-linear -searches `Emcee` and `Zeus`. Like PySwarms, these also need a good starting point, and are generally less effective at +searches `Emcee` and `Zeus`. These also need a good starting point, and are generally less effective at lens modeling than Nautilus. I've included an example runs of Emcee and Zeus below, where the model is set up using `UniformPriors` to give