diff --git a/config/non_linear/mle.yaml b/config/non_linear/mle.yaml index 365b911f..2cc331c7 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 1fbe88cb..079c0086 100644 --- a/config/non_linear/nest.yaml +++ b/config/non_linear/nest.yaml @@ -2,7 +2,7 @@ # - 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. @@ -90,48 +90,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. - diff --git a/scripts/guides/modeling/searches.py b/scripts/guides/modeling/searches.py index 1ea71a97..95e94867 100644 --- a/scripts/guides/modeling/searches.py +++ b/scripts/guides/modeling/searches.py @@ -26,7 +26,7 @@ **Dynesty**: A nested sampling algorithm that is effective for modeling, with a lot of customization. **Emcee**: An ensemble MCMC sampler that is commonly used in Astronomy and Astrophysics. **Zeus**: An ensemble MCMC slice sampler that is the most effective MCMC method for modeling. -**PySwarms**: A particle swarm optimization (PSO) algorithm that is a maximum likelihood estimator (MLE) method. +**LBFGS**: A quasi-Newton optimization algorithm that is a maximum likelihood estimator (MLE) method. **Start Point**: An API that allows the user to specify the start-point of a model-fit, which is useful for MCMC and MLE methods. **Search Cookbook**: A cookbook that documents all searches available in **PyAutoFit**, including those not documented here. @@ -41,7 +41,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 autogalaxy as ag @@ -172,34 +171,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 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 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 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 model very fast. - -In our experience, the parameter spaces fitted by 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 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, ) """ @@ -296,8 +282,8 @@ """ __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. +There are a number of other searches supported by **PyAutoFit** and therefore which can be used, which are not +explictly documented here. 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 f6365605..8f6826d6 100644 --- a/scripts/guides/plot/examples/searches.py +++ b/scripts/guides/plot/examples/searches.py @@ -14,7 +14,6 @@ - **Setup:** Sets up a dataset and model which we will perform quick model-fits to for illustration. - **Dynesty:**: Plots results of the nested sampling method Dynesty. - **Emcee:**: Plots results of an Emcee fit (e.g. cornerplot). -- **PySwarms:**: Plots results of a PySwarms fit (e.g. contour). - **Zeus:**: Plots results of a Zeus fit (e.g. cornerplot). - **GetDist:**: Plot results of any MCMC / nested sampler non-linear search using the library GetDist. @@ -454,39 +453,6 @@ plt.show() -""" -__PySwarms__ - -For PySwarms MLE searches, the `search_internal` attribute contains the pyswarms optimizer object -which can be used for search-specific visualization. -""" -search_internal = result.search_internal - -if search_internal is not None: - 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 Exception: - pass - """ __GetDist__ diff --git a/scripts/howtogalaxy/chapter_optional/tutorial_searches.py b/scripts/howtogalaxy/chapter_optional/tutorial_searches.py index 59ba0154..786b59b6 100644 --- a/scripts/howtogalaxy/chapter_optional/tutorial_searches.py +++ b/scripts/howtogalaxy/chapter_optional/tutorial_searches.py @@ -14,7 +14,7 @@ __Contents__ **Nested Sampling:** Customize Nautilus settings and use Dynesty as an alternative nested sampler. -**Optimizers:** Use maximum likelihood optimizers like PySwarms for fast but less robust fitting. +**Optimizers:** Use maximum likelihood optimizers like LBFGS for fast but less robust fitting. **MCMC:** Use Markov Chain Monte Carlo methods like Emcee for parameter estimation. """ @@ -202,89 +202,19 @@ error estimates (perhaps this is our final 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 galaxy 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 model and therefore needs a the initial swarm of -particles to surround the highest likelihood 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 simulation 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 modeling search chaining (see HowToGalaxy chapter 3). However, even in such -circumstances, we have found that is often unrealible and often infers a local maxima. -""" -lens_bulge = af.Model(ag.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(ag.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(ag.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(ag.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(ag.Galaxy, redshift=0.5, mass=mass, shear=shear) -source = af.Model(ag.Galaxy, redshift=1.0, bulge=bulge) - -model = af.Collection(galaxies=af.Collection(lens=lens, source=source)) - -search = af.PySwarmsLocal( - path_prefix=Path("howtogalaxy", "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 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 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 models are often too complex for optimizers without careful +setup of initialization priors. __MCMC__