Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 0 additions & 44 deletions config/non_linear/mle.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
47 changes: 1 addition & 46 deletions config/non_linear/nest.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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.

38 changes: 12 additions & 26 deletions scripts/guides/modeling/searches.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand All @@ -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
Expand Down Expand Up @@ -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,
)

"""
Expand Down Expand Up @@ -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.
Expand Down
34 changes: 0 additions & 34 deletions scripts/guides/plot/examples/searches.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand Down Expand Up @@ -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__

Expand Down
90 changes: 10 additions & 80 deletions scripts/howtogalaxy/chapter_optional/tutorial_searches.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.
"""

Expand Down Expand Up @@ -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__

Expand Down