Optimization and sampling¶
As mentioned in the Walkthrough, PERI works by fitting a
generative model to data. In the Walkthrough, we obtained
an initial guess for the particle positions and radii from centroid methods. We
used this initial guess to create a state st
with all the components of the
model. We now need to optimize the model.
First, the centroid methods do not provide a guess for global parameters such
as the illumination field or the point-spread function. Second, the centroid
methods frequently miss or mis-feature many particles. Third, we would like to
completely optimize the model to get maximally-accurate position and radii
measurements. In this section, we’ll discuss how to optimize a
peri.states.ImageState
and accomplish all three of these goals.
Most of these steps are implemented in the peri.runner
convenience
functions.
Initial Optimization¶
The goal of the initial optimization is to fit the image accurately enough to identify any missing particles. The initial guess for component values such as the illumination and background fields are usually wildly off. Moreover, centroid methods frequently return wildly inaccurate particle positions. Both the global parameters such as the illumination and the local parameters such as particle positions and radii need to be reasonable before we can find the missing particles.
peri.opt.optimize.burn()
effectively runs this initial optimization. To do this
initial optimization, simply run:
import peri.opt.optimize as opt
opt.burn(st, mode='burn', n_loop=4, desc='', fractol=0.1)
This will usually optimize a state sufficiently enough to add missing particles
to the image. If it does not, you can set n_loop
to a larger value such as
6 or 10. Additionally, opt.burn
saves a copy of the state every so often,
by calling peri.states.save(st, desc='burning')
. If you do not want to
save the state every so oftern, set desc=None
in opt.burn
; if you want
to save with a different name then set desc
to whatever you want; it will
be passed through to peri.states.save()
.
Briefly, what does burn()
do? The state are optimized by
essentially curve-fitting with a Levenberg-Marquardt algorithm.
However, a typical state has too many parameters to optimize the entire state
at once. To deal with this, opt.burn
optimizes the parameters in groups, first
optimizing the global parameters such as the illumination and background, then
optimizing the particle parameters. Since these two groups of parameters are
coupled, opt.burn
then alternates back and forth between globals and
particles n_loop
times or until a sufficient convergence is met – here
if the error changes by less than 10% of its best value (fractol=0.1
) [1].
Finally, we’ve found empirically that it’s best to avoid optimizing some
parameters of the model (such as the point-spread function and the particle
zscale) until the end. Calling burn()
with
mode='burn'
combines all these nuances into one convenience function [2].
[1] | If burn doesn’t converge sufficiently and instead stops because
it has iterated for n_loop , then it will log a warning. For this
initial optimization we’ll ignore it, since we’ll keep optimizing the
state later. However, if you get a warning like this in the later steps,
it might be worthwile to heed the warning and re-optimize the state. |
[2] | In addition, when the point-spread function or illumination are far
from their correct values, the particle radii tend to drift to bad
values before slowly coming back to the correct values. If you have a
pretty good initial guess for the radii, you can speed this up by
having opt.burn ignore the particle radii with the flag
include_rad=False . |
If you want more details on how opt.burn
functions, see the documentation
or code. If your specific image or model is not optimized well by opt.burn
,
or you want additional functionality, then you should look at these functions
which opt.burn
calls or uses:
peri.opt.optimize.do_levmarq()
- Levenberg-Marquardt (LM) optimization on whatever parameter groups passed, additionall optimized for large parameter spaces.
peri.opt.optimize.do_levmarq_all_particle_groups()
- LM optimization on all the particles in the image.
peri.opt.optimize.do_levmarq_particles()
- LM optimization on a select number of particles.
peri.opt.optimize.LMGlobals
- The class that
peri.opt.optimize.do_levmarq()
calls to do its optimization. Has more options and attributes which are useful for checking convergence.peri.opt.optimize.LMParticleGroupCollection
- The class that
peri.opt.optimize.do_levmarq_all_particle_groups()
calls to do its optimization. Has more options and attributes which are useful for checking convergence.peri.opt.optimize.LMParticles
- The class that both
peri.opt.optimize.do_levmarq_particles()
andperi.opt.optimize.LMParticleGroupCollection
calls to do their optimization. Has more options and attributes which are useful for checking convergence.peri.opt.optimize.LMAugmentedState
- Like
LMGlobals
but also allows for effective parameters such as an overall radii scale or a radii scale that changes withz
.peri.opt.optimize.LMEngine
- The workhorse optimizer base class, called by
LMGlobals
andLMParticles
Add-subtract¶
After the initial optimization we can add any missing particles and remove any particles that shouldn’t be there. To do this, run:
import peri.opt.addsubtract as addsub
num_changed, removed_positions, added_positions = addsub.add_subtract(st,
rad='calc', min_rad='calc', max_rad='calc', invert=True,
max_npart='calc')
This function adds missing particles to the image by centroid-featuring the
residuals, with invert
the same as for the initial centroid featuring –
set invert=True
if the particles are dark on a bright background; False
otherwise. In the residuals image, missing particles stick out like sore thumbs
and are easy to find. The function adds a particle at this position with radius
rad
; setting rad='calc'
makes the function choose the radius internally
as the median radius of all the current particles.
More commonly however, two particles are initially featured as one. The initial
optimization will then split the difference by placing this one particle at a
position between the two particles and giving it a large radius. As a result,
the group of particles gets missed by the centroid featuring and particles are
not added. To combat this, the add_subtract()
removes particles that
have a suspiciously large or small radii values, as determined by min_rad
and max_rad
. (Setting these two to 'calc'
uses the cutoffs at the
median radius +/- 15 radii standard deviations.) With the incorrect large
particles removed, the missing particles can be featured. The function
repeatedly removes bad particles and adds missing particles until either no
change is made or it has iterated over the maximum number of loops.
Main Optimization¶
After adding all the particles, it’s time to completely optimize the state. In my experience, sometimes adding particles causes the globals and the old particle positions to no longer be correct. To deal with this, run something like
opt.burn(st, mode='burn', n_loop=6, desc='', fractol=1e-2)
This usually sets the illumination and particle positions to reasonable values.
At this point, it’s time to optimize all the state including the point-spread
function, which we have so far ignored. We can include this with rest of the
parameters with burn()
again:
opt.burn(st, mode='polish', n_loop=6, desc='')
What does this do? First, especially if the initial guess for the point-spread
function was correct, running another optimization with mode='burn'
keeps
the point-spread function from drifting to a bad space because of its strong
coupling with the illumination field. Setting mode='polish'
then causes
burn to optimize everything, alternating between an iteration of optimizing all
the global parameters (including the PSF) and an iteration of optimizing all
the particle positions. Similar to mode='burn'
, setting mode='polish'
saves the state after each iteration by calling
peri.states.save(st, desc='polishing')
; you can set desc
to something
else if you’d like.
Achieving the best-possible state¶
Sometimes, after all this, particles are still missing or the fit is still not
perfect. There are still a few more tricks in the peri
package to fix these
problems. These tricks are incorporated in the peri.runner.finish_state()
for convenience. In case you find a better protocol for your images, here are
the tricks below.
Adding tough missing particles¶
Sometimes one pass of add_subtract()
is not enough
to find all the missing particles, or running the secondary optimizations
reveals that more particles are missing. In these cases, running another
add_subtract()
usually fixes the problem and gets
all the particles. However, sometimes there are particles that the normal
add_subtract()
just can’t seem to get right, such
as particles on the edge or tough clusters of particles. For these cases, there
is another function in the peri.opt.addsubtract
module:
num_added, added_positions = addsub.add_subtract_locally(st)
Briefly, add_subtract_locally()
looks for poorly-
fit regions where the residuals deviate from white Gaussian noise, with the
size of the region roughly set by the optional parameter filter_size
, and
the threshold for badness set by sigma_cutoff
. The function then removes
all the particles in that region and re-adds them based on centroid featuring
again. Since add_subtract_locally()
removes all the
particles in the region, it’s best not to use it until the image is fairly well
fit. Otherwise, the function will attempt to remove nearly all the particles in
the image and re-add them, which takes a long time and will probably fail. That
being said, this function is excellent at fixing doubly-featured particle and
at identifying particles at the edge of or slightly outside of the image. You
can improve its chances of identifying particles at the very edge of the image
by passing a minmass
parameter; I find that minmass=0
frequently works
for particles at the edge of a well-fit state.
Additional Optimizations¶
Occasionally the number of optimization loops isn’t enough to completely
optimize a state. You can try to fix this by running a few more loops of
burn()
with mode='burn'
or mode='polish'
,
depending on whether the illumination is sufficiently far from the minimum as
to bias the PSF. But especially for large states, sometimes this is not enough.
In this case use
opt.finish(st)
burn()
makes some numerical approximations in the
least-squares fit to accelerate convegence. However, when the fit is near the
global minimum, sometimes burn
can get stuck.
finish()
gets around this, taking a slightly slower
but surer descent to the minimum [3].
[3] | For the algorithm-savvy: A Levenberg-Marquardt algorithm works by
evaluating the derivative of the each residuals pixel with respect to
each parameter. A typical microscope image with peri has a few
million pixels and a few thousand parameters, making these derivatives
much too big to store in memory. burn gets around this by treating
the particles separately from the microscope parameters, and only
fitting the microscope parameters with a random subset of the image. But
taking this random subset causes burn to get slightly stuck with the
microscope parameters. finish gets around this by using the entire
image and alternating between fitting only a small subset of the
parameters at once. In addition, by default finish treats the PSF
separately, since it is more difficult to fit. |
Since I like to be sure that I’m at a global minimum, I always run a few extra
loops of opt.burn
with mode='polish'
or of opt.finish
no matter
what.
What if the optimizer gets stuck? If the optimizer is stuck, and you know you are not at the minimum, then you can individually optimize certain parameters. For instance, if you know the PSF is not correct based on the way the residuals looks, you can specifically optimize the PSF by doing this:
opt.do_levmarq(st, st.get('psf').params)
or whatever global component you think is poorly optimized. If the error of the
state decreases significantly, then the state was not at the global minimum and
should be sent through another few loops of opt.burn
or opt.finish
.
When is the state optimized?¶
PERI relies on finding a global minimum of the fit. If the fit is not correct, then obviously your extracted parameters such as positions and radii will not be correct. How can you check if the state is optimized? Below are a few things we check to see if a state is optimized. You can find many of these detailed in the Supplemental Information for our paper.
Checking optimization with the OrthoManipulator¶
The best tool for checking optimization is the peri.viz.interaction.OrthoManipulator
:
from peri.viz.interaction import OrthoManipulator
OrthoManipulator(st)
This will pull up an interactive view of the state, with the data in one panel
and the state model in another. Pressing Q
will cycle the view in the
second panel through the reconstructed model, the fit residuals, and individual
model components. To see if the fit is good, look at the fit residuals. Are
there missing particles, both in the middle of the frame and near the edges?
Can you see shadows of particles? If so, then the state is not optimized. In
contrast, if the residuals are nearly perfect white Gaussian noise, then you’re
done.
The OrthoManipulator
has a lot of additional
functionality, including a view of the Fourier transform of the residuals and
the ability to add, remove, or optimize individual particles interactively. Try
it!
Checking optimization by running more optimization¶
Another way to check is simply to run more loops of opt.burn
or
opt.finish
. If the error or the parameters you care about change
significantly, then you probably needed to run more optimization loops. If not,
then you were near the minimum. While doing this for every image is probably
impractical, you can check a few images or a smaller section of an image to
see if your protocol is good.
Seeing if the fitted values are reasonable¶
Frequently it’s possible to tell if the fit is good simply by looking at the
parameters themselves. Do the particle radii vary systematically with x
,
y
, or z
? If so, then the image is probbaly not at a good fit. We’ve
found that variations in x
or y
tend to be due to imperfections in the
ILM, which varies strongly in these directions for us, and variations in z
tend to be due to imperfections in the PSF, due to the increased aberration
with depth. Note that this might not just be a case of a poor fit – the
complexity of the model could be insufficient. You might need to use a more
realistic PSF or use a higher order for the ILM.
You can do similar checks by looking at either the fitted parameters of the PSF
and other components, or the actual fields themselves using the
OrthoManipulator
or OrthoViewer
.
As an aside, we don’t find it terribly useful to check if the residuals are at
the expected level of the noise. If you somehow knew exactly what the noise
level was, then you could check that st.residuals.std()
is what it should
be. However, the difference between a good fit and a poor fit can be one-tenth
of a percent (i.e. 1e-3) of the residuals. It is highly unlikely that you know
the level of the noise to that precision – the noise level can vary by more
than that from frame-to-frame in a movie due to photobleaching or laser power
fluctuations.
Comparing across Images¶
Finally, you can compare parameters across images. If you featured multiple
images the same way, and the global parameters differ considerably (by
considerably more than the Cramer-Rao Bound), then the state is either not
fully optimized or the model is incomplete. The same applies if the particle
radii fluctuate considerably from frame-to-frame. You can check this easily
with peri.test.track.calculate_state_radii_fluctuations()
.
Speeding this process up¶
Doing this process from start to finish can take a considerable amount of time. In addition to the parallelization methods mentioned in Parallel</parallel>, here are a collection of several tricks to finding a good fit faster.
Using a Good Initial Guess; The Runner Functions¶
The best method for speeding up the featuring is to use a good initial guess. If you know the ILM or the PSF accurately, then you can certainly save time by avoiding the initial fits in Initial Optimization, and possibly even the final fits in Main Optimization.
peri
has convenience functions to use previously-optimized global paramters
in fitting an image. These (and others) are located in the runner
module.
For instance, if you have a previously featured state saved as state_name
,
this will feature a new image '1.tif'
:
from peri import runner
feature_diam = 5 #or whatever feature_diam is best for centroid methods
actual_rad = 5.38 #the actual radius in pixels of the particles
st = runner.get_particles_featuring(feature_diam, state_name=state_name,
im_name='1.tif', actual_rad=actual_rad, do_polish=True)
runner.get_particles_featuring
takes all of the global parameters from the
state state_name
, switches the image, and re-features an initial guess with
centroid methods. It then optimizes the particle positions and radii before
returning. Setting do_polish
to True will automatically run an
opt.burn(st, mode='polish')
on both the globals and particles before
returning the state (takes more time), but this can be omitted for speed. The
state is automatically saved at several points. Similar functionality is
provided by some other of the runner
functions – for instance,
runner.translate_featuring
if the particle positions haven’t moved much
between the loaded state and the new image.
Fitting a small image¶
The larger the image is, the longer it takes to fit. Fitting a small image
considerably speeds up the fit. You can change the region of the fit by setting
the peri.util.Tile
of the image, as described in the
Walkthrough.
Fitting a small image is useful to get a good estimate of global parameters,
especially the point-spread function. Since the exact point-spread functions
included in peri change only with z
, fitting a small portion of the image
in x
and y
but over the full z
extent will still give an accurate
PSF.
We highly encourage you do fit a small image very well to find a good PSF. The
PSF is difficult to optimize (its optimization space is far from a simple
quadratic, and there are slow directions in the fit). Highly-optimizing a small
state to get an accurate PSF will do more than save a lot of time later. For
larger states the optimizer can even get stuck and terminate, thinking it is at
a good fit when it reality the PSF is far from the minimum, which can severely
bias your fits. Make a small image and optimize it overnight – say, 50-100
loops of opt.burn
with mode='polish'
to ensure that you’ve found the
global minimum. You might even want to alternate a loop of burn with a direct
minimization of the PSF, like so:
import numpy as np
state_vals = [] #storing to check at the end
for i in xrange(50): #or another big number
opt.burn(st, mode='polish', n_loop=1)
opt.do_levmarq(st, st.get('psf').params)
state_vals.append(np.copy(st.state[st.params]))
When it finishes, check that the parameters have stopped changing by plotting
them. For instance, to check the parameter psf-alpha
:
import matplotlib.pyplot as plt
index = st.params.index('psf-alpha')
plt.plot(state_vals[:,index])
You should see it smoothly approach a constant value. If it doesn’t look converged, then keep optimizing. If you change your imaging conditions – the index of refraction of the solvent or the microscope and lens – then you will need to do this again.
Sacrificing Precision for Speed¶
peri
is designed to extract information at the maximum possible accuracy.
It does this by finding the best fit of an accurate model to the data. If you
don’t need the maximal possible accuracy, then running peri
to completion
is overkill. For instance, if you just want to distinguish which size a
particle is in a bidisperse suspension, finding the particle radii accurately
to 1 nm is not necessary.
If this is the case, you can save some time by running less optimization loops
or not worrying about finding every last particle. You might also be able to
save time by using a less accurate model – for instance, you could use an ILM
of lower order to create less parameters to fit, or a less accurate PSF to
decrease the execution time for one model generation. You can find some of
these inexact PSFs in peri.comp.psfs
, along with a description of how well
they work in the paper’s Supplemental Information.
Inventing a new algorithm for fitting in high-dimensional spaces¶
Please do this.