peri.opt

peri.opt.optimize.burn

peri.opt.optimize.burn(s, n_loop=6, collect_stats=False, desc=”, rz_order=0, fractol=0.0001, errtol=0.01, mode=’burn’, max_mem=1000000000.0, include_rad=True, do_line_min=’default’, partial_log=False, dowarn=True)

Optimizes all the parameters of a state.

Burns a state through calling LMParticleGroupCollection and LMGlobals/ LMAugmentedState.

Parameters:
  • s (peri.states.ImageState) – The state to optimize
  • n_loop (Int, optional) – The number of times to loop over in the optimizer. Default is 6.
  • collect_stats (Bool, optional) – Whether or not to collect information on the optimizer’s performance. Default is False.
  • desc (string, optional) – Description to append to the states.save() call every loop. Set to None to avoid saving. Default is ”, which selects one of ‘burning’, ‘polishing’, ‘doing_positions’
  • rz_order (Int, optional) – Set to an int > 0 to optimize with an augmented state (R(z) as a global parameter) vs. with the normal global parameters; rz_order is the order of the polynomial approximate for R(z). Default is 0 (no augmented state).
  • fractol (Float, optional) – Fractional change in error at which to terminate. Default 1e-4
  • errtol (Float, optional) – Absolute change in error at which to terminate. Default 1e-2
  • mode ({'burn', 'do-particles', or 'polish'}, optional) – What mode to optimize with. * ‘burn’ : Your state is far from the minimum. * ‘do-particles’ : Positions far from minimum, globals well-fit. * ‘polish’ : The state is close to the minimum. ‘burn’ is the default. Only polish will get to the global minimum.
  • max_mem (Numeric, optional) – The maximum amount of memory allowed for the optimizers’ J’s, for both particles & globals. Default is 1e9, i.e. 1GB per optimizer.
  • do_line_min (Bool or 'default', optional) – Set to True to do an additional, third optimization per loop which optimizes along the subspace spanned by the last 3 steps of the burn()’s trajectory. In principle this should signifi- cantly speed up the convergence; in practice it sometimes does, sometimes doesn’t. Default is ‘default’, which picks by mode: * ‘burn’ : False * ‘do-particles’ : False * ‘polish’ : True
  • dowarn (Bool, optional) – Whether to log a warning if termination results from finishing loops rather than from convergence. Default is True.
Returns:

Dictionary of convergence information. Contains whether the optimization has converged (key 'converged'), the values of the state after each loop (key 'all_loop_values'). The values of the state’s parameters after each part of the loop: globals, particles, linemin. If collect_stats is set, then also contains lists of termination dicts from globals, particles, and line minimization (keys 'global_stats', 'particle_stats', and 'line_stats’).
Return type:

dictionary

Notes

Proceeds by alternating between one Levenberg-Marquardt step optimizing the globals, one optimizing the particles, and repeating until termination.

In addition, if do_line_min is True, at the end of each loop step an additional optimization is tried along the subspaced spanned by the steps taken during the last 3 loops. Ideally, this changes the convergence from linear to quadratic, but it doesn’t always do much.

Each of the 3 options proceed by optimizing as follows: * burn : lm.do_run_2(), lp.do_run_2(). No psf, 2 loops on lm. * do-particles : lp.do_run_2(), scales for ilm, bkg’s * polish : lm.do_run_2(), lp.do_run_2(). Everything, 1 loop each. where lm is a globals LMGlobals instance, and lp a LMParticleGroupCollection instance.

peri.opt.optimize.finish

peri.opt.optimize.finish(s, desc=’finish’, n_loop=4, max_mem=1000000000.0, separate_psf=True, fractol=1e-07, errtol=0.001, dowarn=True)

Crawls slowly to the minimum-cost state.

Blocks the global parameters into small enough sections such that each can be optimized separately while including all the pixels (i.e. no decimation). Optimizes the globals, then the psf separately if desired, then particles, then a line minimization along the step direction to speed up convergence.

Parameters:
  • s (peri.states.ImageState) – The state to optimize
  • desc (string, optional) – Description to append to the states.save() call every loop. Set to None to avoid saving. Default is ‘finish’.
  • n_loop (Int, optional) – The number of times to loop over in the optimizer. Default is 4.
  • max_mem (Numeric, optional) – The maximum amount of memory allowed for the optimizers’ J’s, for both particles & globals. Default is 1e9.
  • separate_psf (Bool, optional) – If True, does the psf optimization separately from the rest of the globals, since the psf has a more tortuous fit landscape. Default is True.
  • fractol (Float, optional) – Fractional change in error at which to terminate. Default 1e-4
  • errtol (Float, optional) – Absolute change in error at which to terminate. Default 1e-2
  • dowarn (Bool, optional) – Whether to log a warning if termination results from finishing loops rather than from convergence. Default is True.
Returns:

Information about the optimization. Has two keys: 'converged', a Bool which of whether optimization stopped due to convergence (True) or due to max number of iterations (False), and 'loop_values', a [n_loop+1, N] numpy.ndarray of the state’s values, at the start of optimization and at the end of each loop, before the line minimization.
Return type:

dictionary

peri.opt.optimize

peri.opt.optimize functions

peri.opt.optimize.do_levmarq(s, param_names, damping=0.1, decrease_damp_factor=10.0, run_length=6, eig_update=True, collect_stats=False, rz_order=0, run_type=2, **kwargs)

Runs Levenberg-Marquardt optimization on a state.

Convenience wrapper for LMGlobals. Same keyword args, but the defaults have been set to useful values for optimizing globals. See LMGlobals and LMEngine for documentation.

See also

do_levmarq_particles()
Levenberg-Marquardt optimization of a specified set of particles.
do_levmarq_all_particle_groups()
Levenberg-Marquardt optimization of all the particles in the state.
LMGlobals()
Optimizer object; the workhorse of do_levmarq.
LMEngine()
Engine superclass for all the optimizers.
peri.opt.optimize.do_levmarq_all_particle_groups(s, region_size=40, max_iter=2, damping=1.0, decrease_damp_factor=10.0, run_length=4, collect_stats=False, **kwargs)

Levenberg-Marquardt optimization for every particle in the state.

Convenience wrapper for LMParticleGroupCollection. Same keyword args, but I’ve set the defaults to what I’ve found to be useful values for optimizing particles. See LMParticleGroupCollection for documentation.

See also

do_levmarq_particles()
Levenberg-Marquardt optimization of a specified set of particles.
do_levmarq()
Levenberg-Marquardt optimization of the entire state; useful for optimizing global parameters.
LMParticleGroupCollection()
The workhorse of do_levmarq.
LMEngine()
Engine superclass for all the optimizers.
peri.opt.optimize.do_levmarq_particles(s, particles, damping=1.0, decrease_damp_factor=10.0, run_length=4, collect_stats=False, max_iter=2, **kwargs)

Levenberg-Marquardt optimization on a set of particles.

Convenience wrapper for LMParticles. Same keyword args, but the defaults have been set to useful values for optimizing particles. See LMParticles and LMEngine for documentation.

See also

do_levmarq_all_particle_groups()
Levenberg-Marquardt optimization of all the particles in the state.
do_levmarq()
Levenberg-Marquardt optimization of the entire state; useful for optimizing global parameters.
LMParticles()
Optimizer object; the workhorse of do_levmarq_particles.
LMEngine()
Engine superclass for all the optimizers.

peri.opt.optimize classes

class peri.opt.optimize.LMEngine(damping=1.0, increase_damp_factor=3.0, decrease_damp_factor=8.0, min_eigval=1e-13, marquardt_damping=False, transtrum_damping=None, use_accel=False, max_accel_correction=1.0, paramtol=1e-06, errtol=1e-05, exptol=0.001, fractol=1e-06, costol=None, max_iter=5, run_length=5, update_J_frequency=1, broyden_update=True, eig_update=False, eig_update_frequency=3, num_eig_dirs=8, eig_dl=1e-05, broyden_update_frequency=1)

Levenberg-Marquardt engine with all the options from the M. Transtrum J. Sethna 2012 ArXiV paper [1].

Parameters:
  • damping (Float, optional) – The initial damping factor for Levenberg-Marquardt. Adjusted internally. Default is 1.
  • increase_damp_factor (Float, optional) – The amount to increase damping by when an attempted step has failed. Default is 3.
  • decrease_damp_factor (Float, optional) – The amount to decrease damping by after a successful step. Default is 8. increase_damp_factor and decrease_damp_factor must not have all the same factors.
  • min_eigval (Float scalar, optional, <<1.) – The minimum eigenvalue to use in inverting the JTJ matrix, to avoid degeneracies in the parameter space (i.e. ‘rcond’ in np.linalg.lstsq). Default is 1e-12.
  • marquardt_damping (Bool, optional) – Set to False to use Levenberg damping (damping matrix proportional to the identiy) instead of Marquardt damping (damping matrix proportional to the diagonal terms of JTJ). Default is False.
  • transtrum_damping (Float or None, optional) – If not None, then clips the Marquardt damping diagonal entries to be at least transtrum_damping. Default is None.
  • use_accel (Bool, optional) – Set to True to incorporate the geodesic acceleration term from M. Transtrum J. Sethna 2012. Default is False.
  • max_accel_correction (Float, optional) – Acceleration corrections bigger than max_accel_correction times the normal LM step are viewed as bad steps, causing a decrease in damping. Default is 1.0. Only applies to the do_run_1 method.
  • paramtol (Float, optional) – Algorithm has converged when the none of the parameters have changed by more than paramtol. Default is 1e-6.
  • errtol (Float, optional) – Algorithm has converged when the error has changed by less than errtol after 1 step. Default is 1e-6.
  • exptol (Float, optional) – Algorithm has converged when the expected change in error is less than exptol. Default is 1e-3.
  • fractol (Float, optional) – Algorithm has converged when the error has changed by a fractional amount less than fractol after 1 step. Default is 1e-6.
  • costol (Float, optional) – Algorithm has converged when the cosine of the angle between (residuals projected onto the model manifold) and (the residuals) is < costol. Default is None, i.e. doesn’t check the cosine (since it takes a bit of time).
  • max_iter (Int, optional) – The maximum number of iterations before the algorithm stops iterating. Default is 5.
  • update_J_frequency (Int, optional) – The frequency to re-calculate the full Jacobian matrix. Default is 2, i.e. every other run.
  • broyden_update (Bool, optional) – Set to True to do a Broyden partial update on J after each step, updating the projection of J along the parameter change direction. Cheap in time cost, but not always accurate. Default is False.
  • eig_update (Bool, optional) – Set to True to update the projection of J along the most stiff eigendirections of JTJ. Slower than broyden but more accurate & useful. Default is False.
  • num_eig_dirs (Int, optional) – If eig_update == True, the number of eigendirections to update when doing the eigen update. Default is 4.
  • eig_update_frequency (Int, optional) – If eig_update, the frequency to do this partial update. Default is 3.
  • broyden_update_frequency (Int, optional) – If broyden_update, the frequency to do this partial update. Default is 1.
J

numpy.ndarray – The calculated J for Levenberg-Marquardt. Starts as None

JTJ

numpy.ndarray – The approximation to the fit Hessian; np.dot(J, J.T)

damping

numpy.ndarray – The current damping vector for the parameters.

_num_iter

Int – The number of iterations ran in the current cycle. Don’t touch

reset(new_damping=None)

Resets all counters etc so a new run can commence.

do_run_1()

Method 1 for optimization.

do_run_2()

Method 2 for optimization

do_internal_run()

Additional, slight optimization once J has been calculated

find_LM_updates(grad, do_correct_damping=True, subblock=None)

Returns the Levenberg-Marquardt step.

increase_damping()

Increases damping

decrease_damping(undo_decrease=False)

Decreases damping or undoes a previous decrease.

update_param_vals(new_vals, incremental=False)

Updates the current set of parameter values and previous values, sets a flag to re-calculate J.

calc_model_cosine(decimate=None)

Calculates the cosine of the residuals with the current model tangent space J

get_termination_stats(get_cos=True)

Returns a dict of termination statistics

check_completion()

Checks if the algorithm has found a satisfactory minimum

check_termination()

Checks if the algorithm should terminate

update_J()

Updates J, JTJ

calc_grad()

Calculates the gradient of the cost w.r.t. the parameters.

update_Broyden_J()

Execute a Broyden update of J

update_eig_J()

Execute an eigen update of J

calc_accel_correction(damped_JTJ, delta0)

Calculates the geodesic acceleration correction to the standard Levenberg-Marquardt step.

update_select_J(blk)

Updates J only for certain parameters, described by the boolean mask blk.

See also

LMFunction
Levenberg-Marquardt (LM) optimization for a user- supplied function.
LMGlobals
LM optimization designed for optimizing state globals.
LMParticles
LM optimization designed for optimizing state particles.
LMParticleGroupCollection
LM optimization on all particles in a state.
LMAugmentedState
LMGlobals with additional R(z) parameters.
LMOptObj
…

Notes

There are 3 different options for optimizing:

  • do_run_1():
    Checks to calculate full, Broyden, and eigen J, then tries a step. If the step is accepted, decreases damping; if not, increases.
  • do_run_2():
    Checks to calculate full, Broyden, and eigen J, then tries a step with the current damping and with a decreased damping, accepting whichever is lower. Decreases damping iff the lower damping is better. It then calls do_internal_run() (see below). Rejected steps result in increased damping until a step is accepted. Checks for full, Broyden, and eigen J updates.
  • do_internal_run():
    Checks for Broyden and eigen J updates only, then uses pre-calculated J, JTJ, etc to evaluate LM steps. Does not change damping during the run. Does not check do update the full J, but does check for Broyden, eigen updates. Does not work if J has not been evaluated yet.

Whether to update the full J is controlled by update_J_frequency only, which only counts iterations of do_run_1() and do_run_2(). Partial updates are controlled by *_update_frequency, which counts internal runs in do_internal_run and full runs in do_run_1.

So, if you want a partial update every other run, full J the remaining, this would be:

  • do_run_1(): update_J_frequency=2, partial_update_frequency=1
  • do_run_2(): update_J_frequency=1, partial_update_frequency=1, run_length=2

Partial Updates:

  • Broyden update : an update to J (and then JTJ) by approximating the
    derivative of the model as the finite difference of the last step. (rank-1)
  • Eigen update : a rank-num_eig_dirs update to J (and then JTJ) by
    finite-differencing with eig_dl along the highest num_eig_dirs eigendirections.

Damping:

  • marquardt : Damp proportional to the diagonal elements of JTJ
  • transtrum : Marquardt damping, clipped to be at least a certain number
  • default : (levenberg) Damp using something proportional to the identity

References

[1]M. Transtrum and J. Sethna, “Improvements to the Levenberg- Marquardt algorithm for nonlinear least-squares minimization,” ArXiV preprint arXiv:1201.5885 (2012)
calc_J()

Updates self.J, returns nothing

calc_accel_correction(damped_JTJ, delta0)

Geodesic acceleration correction to the LM step.

Parameters:
  • damped_JTJ (numpy.ndarray) – The damped JTJ used to calculate the initial step.
  • delta0 (numpy.ndarray) – The initial LM step.
Returns:

corr – The correction to the original LM step.
Return type:

numpy.ndarray
calc_grad()

The gradient of the cost w.r.t. the parameters.

calc_model_cosine(decimate=None, mode=’err’)

Calculates the cosine of the residuals with the model.

Parameters:
  • decimate (Int or None, optional) – Decimate the residuals by decimate pixels. If None, no decimation is used. Valid only with mode=’svd’. Default is None
  • mode ({'svd', 'err'}) – Which mode to use; see Notes section. Default is ‘err’.
Returns:

abs_cos – The absolute value of the model cosine.
Return type:

numpy.float64

Notes

The model cosine is defined in terms of the geometric view of curve-fitting, as a model manifold embedded in a high-dimensional space. The model cosine is the cosine of the residuals vector with its projection on the tangent space: \(cos(phi) = |P^T r|/|r|\) where \(P^T\) is the projection operator onto the model manifold and \(r\) the residuals. This can be calculated two ways: By calculating the projection operator P directly with SVD (mode=`svd`), or by using the expected error if the model were linear to calculate a model sine first (mode=`err`). Since the SVD of a large matrix is slow, mode=`err` is faster.

decimate allows for every nth pixel only to be counted in the SVD matrix of J for speed. While this is n x faster, it is considerably less accurate, so the default is no decimation.

calc_residuals()

returns residuals = data - model.

check_completion()

Returns a Bool of whether the algorithm has found a satisfactory minimum

check_terminate()

Returns a Bool of whether to terminate.

Checks whether a satisfactory minimum has been found or whether too many iterations have occurred.

check_update_J()

Checks if the full J should be updated.

Right now, just updates after update_J_frequency loops

do_internal_run(initial_count=0, subblock=None, update_derr=True)

Takes more steps without calculating J again.

Given a fixed damping, J, JTJ, iterates calculating steps, with optional Broyden or eigendirection updates. Iterates either until a bad step is taken or for self.run_length times. Called internally by do_run_2() but is also useful on its own.

Parameters:
  • initial_count (Int, optional) – The initial count of the run. Default is 0. Increasing from 0 effectively temporarily decreases run_length.
  • subblock (None or np.ndarray of bools, optional) – If not None, a boolean mask which determines which sub- block of parameters to run over. Default is None, i.e. all the parameters.
  • update_derr (Bool, optional) – Set to False to not update the variable that determines delta_err, preventing premature termination through errtol.

Notes

It might be good to do something similar to update_derr with the parameter values, but this is trickier because of Broyden updates and _fresh_J.

do_run_1()

LM run, evaluating 1 step at a time.

Broyden or eigendirection updates replace full-J updates until a full-J update occurs. Does not run with the calculated J (no internal run).

do_run_2()

LM run evaluating 2 steps (damped and not) and choosing the best.

After finding the best of 2 steps, runs with that damping + Broyden or eigendirection updates, until deciding to do a full-J update. Only changes damping after full-J updates.

find_LM_updates(grad, do_correct_damping=True, subblock=None)

Calculates LM updates, with or without the acceleration correction.

Parameters:
  • grad (numpy.ndarray) – The gradient of the model cost.
  • do_correct_damping (Bool, optional) – If self.use_accel, then set to True to correct damping if the acceleration correction is too big. Default is True Does nothing is self.use_accel is False
  • subblock (slice, numpy.ndarray, or None, optional) – Set to a slice or a valide numpy.ndarray to use only a certain subset of the parameters. Default is None, i.e. use all the parameters.
Returns:

delta – The Levenberg-Marquadt step, relative to the old parameters. Size is always self.param_vals.size.
Return type:

numpy.ndarray
find_expected_error(delta_params=’calc’)

Returns the error expected after an update if the model were linear.

Parameters:delta_params ({numpy.ndarray, 'calc', or 'perfect'}, optional) – The relative change in parameters. If ‘calc’, uses update calculated from the current damping, J, etc; if ‘perfect’, uses the update calculated with zero damping.
Returns:The expected error after the update with delta_params
Return type:numpy.float64
get_termination_stats(get_cos=True)

Returns a dict of termination statistics

Parameters:get_cos (Bool, optional) – Whether or not to calcualte the cosine of the residuals with the tangent plane of the model using the current J. The calculation may take some time. Default is True
Returns:
Has keys
delta_vals : The last change in parameter values. delta_err : The last change in the error. exp_err : The expected (last) change in the error. frac_err : The fractional change in the error. num_iter : The number of iterations completed. error : The current error.
Return type:dict
reset(new_damping=None)

Keeps all user supplied options the same, but resets counters etc.

update_Broyden_J()

Execute a Broyden update of J

update_J()

Updates J, JTJ, and internal counters.

update_eig_J()

Execute an eigen update of J

update_function(param_vals)

Takes an array param_vals, updates function, returns the new error

update_param_vals(new_vals, incremental=False)

Updates the current set of parameter values and previous values, sets a flag to re-calculate J.

Parameters:
  • new_vals (numpy.ndarray) – The new values to update to
  • incremental (Bool, optional) – Set to True to make it an incremental update relative to the old parameters. Default is False
update_select_J(blk)

Updates J only for certain parameters, described by the boolean mask blk.

class peri.opt.optimize.LMGlobals(state, param_names, max_mem=1000000000.0, opt_kwargs={}, **kwargs)

Levenberg-Marquardt, optimized for state globals.

Contains alll the options from the M. Transtrum J. Sethna 2012 ArXiV paper. See LMEngine for further documentation.

Parameters:
  • state (peri.states.State) – The state to optimize. Stored as self.state.
  • param_names (List) – List of the parameter names (strings) to optimize over. Stored as self.param_names.
  • max_mem (Numeric, optional) – The maximum memory to use for the optimization; controls pixel decimation. Default is 1e9. Stored as self.max_mem
  • opt_kwargs (Dict, optional) – Dict of **kwargs for opt implementation. Right now only for get_num_px_jtj, i.e. keys of ‘decimate’, ‘min_redundant’. Default is {}. Stored as self.opt_kwargs
num_pix

Int – The number of pixels of the residuals used to calculate J.

set_params(new_param_names, new_damping=None)

Change the parameter names to optimize.

reset(new_damping=None)

Resets counters etc to zero, allowing more runs to commence.

See also

LMEngine
Parent class for all Levenberg-Marquardt (LM) optimization
LMParticles
LM optimization designed for optimizing state particles.
LMParticleGroupCollection
LM optimization on all particles in a state.
LMAugmentedState
LMGlobals with additional R(z) parameters.
do_levmarq
Convenience function for LMGlobals
do_levmarq_particles
Convenience function for optimizing particles
class peri.opt.optimize.LMParticles(state, particles, include_rad=True, **kwargs)

Levenberg-Marquardt, optimized for state globals.

Contains alll the options from the M. Transtrum J. Sethna 2012 ArXiV paper. See LMEngine for further documentation.

Parameters:
  • state (peri.states.ImageState) – The state to optimize
  • particles (numpy.ndarray) – Array of the particle indices to optimize over.
  • include_rad (Bool, optional) – Whether or not to include the particle radii in the optimization. Default is True
param_names

List – The list of the parameter names being optimized.

set_particles(new_particles, new_damping=None)

Change the particle to optimize.

reset(new_damping=None)

Resets counters etc to zero, allowing more runs to commence.

See also

LMEngine
Parent class for all Levenberg-Marquardt (LM) optimization
LMGlobals
LM optimization designed for optimizing state global parameters.
LMParticleGroupCollection
LM optimization on all particles in a state.
LMAugmentedState
LMGlobals with additional R(z) parameters.
do_levmarq
Convenience function for LMGlobals
do_levmarq_particles
Convenience function for optimizing particles

Notes

To prevent the state updates from breaking, this clips the particle rads to [self._MINRAD, self._MAXRAD] and the positions to at least self._MINDIST from the edge of the padded image. These are: * _MINRAD : 1e-3 * _MAXRAD : 2e2 * _MINDIST : 1e-3 For extremely large particles (e.g. larger than _MAXRAD or larger than the pad and barely overlapping the image) these numbers might be insufficient.

class peri.opt.optimize.LMParticleGroupCollection(state, region_size=40, do_calc_size=True, max_mem=1000000000.0, get_cos=False, save_J=False, **kwargs)

Levenberg-Marquardt on all particles in a state.

Convenience wrapper for LMParticles. Generates a separate instance for the particle groups each time and optimizes with that, since storing J for the particles is too large.

Parameters:
  • state (peri.states.ImageState) – The state to optimize
  • region_size (Int or 3-element list-like of ints, optional) – The region size for sub-blocking particles. Default is 40
  • do_calc_size (Bool, optional) – If True, calculates the region size internally based on the maximum allowed memory. Default is True
  • max_mem (Numeric, optional) – The maximum allowed memory for J to occupy. Default is 1e9
  • get_cos (Bool, optional) – Set to True to include the model cosine in the statistics on each individual group’s run, using LMEngine get_termination_stats(). Stored in self.stats. Default is False
  • save_J (Bool) – Set to True to create a series of temp files that save J for each group of particles. Needed for do_internal_run(). Default is False.
Other Parameters:
 
  • Pass any kwargs that would be passed to LMParticles. Stored in
  • self._kwargs for reference.
stats

List – A list of the termination stats for each sub-block of particles

particle_groups

List – A list of particle groups. Element [i] in the list is a numpy.ndarray of the indices in group [i].

region_size

Int or 3-element list-like – The region size of the tiles. If do_calc_size is True, region_size will be the calculated value, which may differ from the input value.

_kwargs

Dict – The **kwargs passed to LMParticles.

reset()

Re-calculate all the groups

do_run_1()

Run do_run_1 for every group of particles

do_run_2()

Run do_run_2 for every group of particles

do_internal_run()

Run do_internal_run for every group of particles

See also

LMParticles
LM optimization designed for optimizing state particles.

Notes

Since storing J for many particles can require a huge amount of memory, this object proceeds by re-initializing a separate LMParticles instance for each group of particles, calculating and then discarding J each time. The calculated J’s can be kept by setting save_J to True, which saves each J for each group in a separate tempfile, located in the current directory. The J’s can then be loaded again to attempt a second step without re-calculating J. However, for a big image this can attempt to store _a_lot_ of temp files, which might be more than the operating system limit (as temp files are always open), which will raise an error. So use with caution. Deleting any references to the temp files by deleting the LMParticleGroupCollection instance will close and remove the temporary files.

class peri.opt.optimize.LMAugmentedState(aug_state, max_mem=1000000000.0, opt_kwargs={}, **kwargs)

Levenberg-Marquardt on an augmented state.

Contains all the options from the M. Transtrum J. Sethna 2012 ArXiV paper. See LMEngine for further documentation.

Parameters:
  • aug_state (:class:peri.optimize.opt.AugmentedState) – The state to optimize. Stored as self.aug_state
  • max_mem (Numeric, optional) – The maximum memory to use for the optimization; controls pixel decimation. Default is 1e9. Stored as self.max_mem.
  • opt_kwargs (Dict, optional) – Dict of **kwargs for opt implementation. Right now only for get_num_px_jtj, i.e. keys of ‘decimate’, min_redundant’. Default is {}. Stored as self.opt_kwargs.
num_pix

Int – The number of pixels of the residuals used to calculate J.

reset(new_damping=None)

Resets the augmented state, counters, etc to zero, allowing more runs to commence.

See also

LMEngine
Parent class for all Levenberg-Marquardt (LM) optimization
LMGlobals
LM optimization on globals without an AugmentedState
AugmentedState
The class used by LMAugmentedState.
LMParticles
LM optimization designed for optimizing state particles.
LMParticleGroupCollection
LM optimization on all particles in a state.
LMAugmentedState
LMGlobals with additional R(z) parameters.
do_levmarq
Convenience function for LMGlobals
do_levmarq_particles
Convenience function for optimizing particles

peri.opt.addsubtract

peri.opt.addsubtract.add_subtract(st, max_iter=7, max_npart=’calc’, max_mem=200000000.0, always_check_remove=False, **kwargs)

Automatically adds and subtracts missing & extra particles.

Operates by removing bad particles then adding missing particles on repeat, until either no particles are added/removed or after max_iter attempts.

Parameters:
  • st (peri.states.State) – The state to add and subtract particles to.
  • max_iter (Int, optional) – The maximum number of add-subtract loops to use. Default is 7. Terminates after either max_iter loops or when nothing has changed.
  • max_npart (Int or 'calc', optional) – The maximum number of particles to add before optimizing the non-psf globals. Default is 'calc', which uses 5% of the initial number of particles.
  • max_mem (Int, optional) – The maximum memory to use for optimization after adding max_npart particles. Default is 2e8.
  • always_check_remove (Bool, optional) – Set to True to always check whether to remove particles. If False, only checks for removal while particles were removed on the previous attempt. Default is False.
Other Parameters:
 
  • invert (Bool, optional) – True if the particles are dark on a bright background, False if they are bright on a dark background. Default is True.
  • min_rad (Float, optional) – Particles with radius below min_rad are automatically deleted. Default is 'calc' = median rad - 25* radius std.
  • max_rad (Float, optional) – Particles with radius above max_rad are automatically deleted. Default is 'calc' = median rad + 15* radius std, but you should change this for your particle sizes.
  • min_edge_dist (Float, optional) – Particles closer to the edge of the padded image than this are automatically deleted. Default is 2.0.
  • check_rad_cutoff (2-element float list.) – Particles with radii < check_rad_cutoff[0] or > check...[1] are checked if they should be deleted (not automatic). Default is [3.5, 15].
  • check_outside_im (Bool, optional) – Set to True to check whether to delete particles whose positions are outside the un-padded image.
  • rad (Float, optional) – The initial radius for added particles; added particles radii are not fit until the end of add_subtract. Default is 'calc', which uses the median radii of active particles.
  • tries (Int, optional) – The number of particles to attempt to remove or add, per iteration. Default is 50.
  • im_change_frac (Float, optional) – How good the change in error needs to be relative to the change in the difference image. Default is 0.2; i.e. if the error does not decrease by 20% of the change in the difference image, do not add the particle.
  • min_derr (Float, optional) – The minimum change in the state’s error to keep a particle in the image. Default is '3sig' which uses 3*st.sigma.
  • do_opt (Bool, optional) – Set to False to avoid optimizing particle positions after adding.
  • minmass (Float, optional) – The minimum mass for a particle to be identified as a feature, as used by trackpy. Defaults to a decent guess.
  • use_tp (Bool, optional) – Set to True to use trackpy to find missing particles inside the image. Not recommended since trackpy deliberately cuts out particles at the edge of the image. Default is False.
Returns:

  • total_changed (Int) – The total number of adds and subtracts done on the data. Not the same as changed_inds.size since the same particle or particle index can be added/subtracted multiple times.
  • added_positions ([N_added,3] numpy.ndarray) – The positions of particles that have been added at any point in the add-subtract cycle.
  • removed_positions ([N_added,3] numpy.ndarray) – The positions of particles that have been removed at any point in the add-subtract cycle.

Notes

Occasionally after the intial featuring a cluster of particles is featured as 1 big particle. To fix these mistakes, it helps to set max_rad to a physical value. This removes the big particle and allows it to be re-featured by (several passes of) the adds.

The added/removed positions returned are whether or not the position has been added or removed ever. It’s possible that a position is added, then removed during a later iteration.

peri.opt.addsubtract.add_subtract_locally(st, region_depth=3, filter_size=5, sigma_cutoff=8, **kwargs)

Automatically adds and subtracts missing particles based on local regions of poor fit.

Calls identify_misfeatured_regions to identify regions, then add_subtract_misfeatured_tile on the tiles in order of size until region_depth tiles have been checked without adding any particles.

Parameters:
  • st (peri.states.State) – The state to add and subtract particles to.
  • region_depth (Int) – The minimum amount of regions to try; the algorithm terminates if region_depth regions have been tried without adding particles.
Other Parameters:
 
  • filter_size (Int, optional) – The size of the filter for calculating the local standard deviation; should approximately be the size of a poorly featured region in each dimension. Best if odd. Default is 5.
  • sigma_cutoff (Float, optional) – The max allowed deviation of the residuals from what is expected, in units of the residuals’ standard deviation. Lower means more sensitive, higher = less sensitive. Default is 8.0, i.e. one pixel out of every 7*10^11 is mis-identified randomly. In practice the noise is not Gaussian so there are still some regions mis- identified as improperly featured.
  • rad (Float or ‘calc’, optional) – The initial radius for added particles; added particles radii are not fit until the end of add_subtract. Default is 'calc', which uses the median radii of active particles.
  • max_iter (Int, optional) – The maximum number of loops for attempted adds at one tile location. Default is 3.
  • invert (Bool, optional) – Whether to invert the image for feature_guess. Default is True, i.e. dark particles on bright background.
  • max_allowed_remove (Int, optional) – The maximum number of particles to remove. If the misfeatured tile contains more than this many particles, raises an error. If it contains more than half as many particles, throws a warning. If more than this many particles are added, they are optimized in blocks of max_allowed_remove. Default is 20.
  • im_change_frac (Float, between 0 and 1.) – If adding or removing a particle decreases the error less than im_change_frac * the change in the image, the particle is deleted. Default is 0.2.
  • min_derr (Float) – The minimum change in the state’s error to keep a particle in the image. Default is '3sig' which uses 3*st.sigma.
  • do_opt (Bool, optional) – Set to False to avoid optimizing particle positions after adding them. Default is True
  • minmass (Float, optional) – The minimum mass for a particle to be identified as a feature, as used by trackpy. Defaults to a decent guess.
  • use_tp (Bool, optional) – Set to True to use trackpy to find missing particles inside the image. Not recommended since trackpy deliberately cuts out particles at the edge of the image. Default is False.
  • max_allowed_remove (Int, optional) – The maximum number of particles to remove. If the misfeatured tile contains more than this many particles, raises an error. If it contains more than half as many particles, throws a warning. If more than this many particles are added, they are optimized in blocks of max_allowed_remove. Default is 20.
Returns:

  • n_added (Int) – The change in the number of particles; i.e the number added - number removed.
  • new_poses (List) – [N,3] element list of the added particle positions.

Notes

Algorithm Description

  1. Identify mis-featured regions by how much the local residuals deviate from the global residuals, as measured by the standard deviation of both.
  2. Loop over each of those regions, and:
    1. Remove every particle in the current region.
    2. Try to add particles in the current region until no more can be added while adequately decreasing the error.
    3. Terminate if at least region_depth regions have been checked without successfully adding a particle.

Because this algorithm is more judicious about chooosing regions to check, and more aggressive about removing particles in those regions, it runs faster and does a better job than the (global) add_subtract. However, this function usually does not work better as an initial add- subtract on an image, since (1) it doesn’t check for removing small/big particles per se, and (2) when the poorly-featured regions of the image are large or when the fit is bad, it will remove essentially all of the particles, taking a long time. As a result, it’s usually best to do a normal add_subtract first and using this function for tough missing or double-featured particles.

peri.opt.addsubtract.feature_guess(st, rad, invert=True, minmass=None, use_tp=False, trim_edge=False, **kwargs)

Makes a guess at particle positions using heuristic centroid methods.

Parameters:
  • st (peri.states.State) – The state to check adding particles to.
  • rad (Float) – The feature size for featuring.
  • invert (Bool, optional) – Whether to invert the image. Default is True, i.e. dark particles
  • minmass (Float or None, optional) – The minimum mass/masscut of a particle. Default is None = calculated internally.
  • use_tp (Bool, optional) – Whether or not to use trackpy. Default is False, since trackpy cuts out particles at the edge.
  • trim_edge (Bool, optional) – Whether to trim particles at the edge pixels of the image. Can be useful for initial featuring but is bad for adding missing particles as they are frequently at the edge. Default is False.
Returns:

  • guess ([N,3] numpy.ndarray) – The featured positions of the particles, sorted in order of decreasing feature mass.
  • npart (Int) – The number of added particles.