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
49 changes: 0 additions & 49 deletions config/non_linear/nest.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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).
9 changes: 0 additions & 9 deletions config/visualize/plots_search.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -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?
Expand Down
36 changes: 11 additions & 25 deletions scripts/guides/modeling/searches.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.

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

"""
Expand Down Expand Up @@ -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.
Expand Down
58 changes: 0 additions & 58 deletions scripts/guides/plot/examples/searches.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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__
Expand Down
92 changes: 11 additions & 81 deletions scripts/howtolens/chapter_optional/tutorial_searches.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.

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