stochrare.rare.ams¶

Rare event algorithms of the Adaptive Multilevel Splitting family¶

There are two kinds of variants of the AMS algorithms. On the one hand, there are “scientific” variants, corresponding to different formulations of the algorithm (e.g. AMS vs TAMS). On the other hand, there are “technical” variants, corresponding to different implementations: for instance, keeping all the trajectories in memory or storing them on disk (necessary for applications to complex systems)

class stochrare.rare.ams.AMS(model, scorefun, initcond=<function AMS.<lambda>>)¶

Bases: object

Original version of the Adaptive Multilevel Splitting Algorithm.

Parameters:

model (stochrare.dynamics.DiffusionProcess1D object (or a subclass of it)) – The dynamical model; so far we are restricted to SDEs of the form \(dX_t = F(X_t, t) + \sqrt{2D}dW_t\). We only use the stochrare.dynamics.DiffusionProcess1D.update() method of the object.
scorefun (function with two arguments) – The score function \(\xi(t, x)\).
initcond (function with no arguments, optional) – Function to generate initial conditions. It can be for instance a constant: lambda: x0, t0 or generate random initial conditions: lambda: np.random.random(), t0

Notes

The algorithm evolves an ensemble of trajectories in an interactive manner, using selection and mutation steps [1] [2] [3] [4]. The algorithm requires two sets \(A\) and \(B\), and a reactive coordinate or score function \(\xi\), measuring the distance between the two. In fact, we require that \(\xi\) vanishes over the boundary of \(A\), and takes unit value over the boundary of \(B\).

Initialization:

The ensemble is initialized by running \(N\) trajectories until they reach set \(A\) or set \(B\).

Selection

Then at each iteration, the maximum value of the score function over each member of the ensemble is computed. The \(q\) trajectories with lowest score function are killed.

Mutation

For each trajectory killed, we pick a random trajectory among the survivors. We clone that trajectory until it reaches the level of the killed trajectory for the first time, then we restart it from that point until it reaches set \(A\) or \(B\).

The algorithm is iterated until all trajectories reach set \(B\).

References

[1]	Cérou and A. Guyader, Stoch. Anal. Appl. 25, 417 (2007)

[2]	Cérou, A. Guyader, T. Lelièvre and D. Pommier J. Chem. Phys. 134, 054108 (2011)

[3]	Rolland, F. Bouchet and E. Simonnet, J. Stat. Phys. 162, 277 (2016)

[4]	C.-E. Bréhier, M. Gazeau, L. Goudenège, T. Lelièvre and M. Rousset, Ann. Appl. Probab. 26, 3559 (2016)

getcrossingtime(level, times, traj)¶

Return the time and position at which the trajectory reaches a given threshold.

Parameters:	level (float) – The threshold. times (numpy.ndarray) – Sampling times for the trajectory. traj (numpy.ndarray) – Position of the system at the sampling times.
Returns:	t, x – The time and position at the crossing point.
Return type:	float, float

getlevel(times, traj)¶

Return the maximum reached by the score function over the trajectory.

Parameters:	times (numpy.ndarray) – Sampling times for the trajectory. traj (numpy.ndarray) – Position of the system at the sampling times.
Returns:	max – The maximum of the score function over the trajectory.
Return type:	float

resample(time, pos, told, xold, **kwargs)¶

Resample a killed trajectory after a given time.

Keyword Arguments:
Parameters:	time (float) – The time from which to resample. pos (float) – The position from which to resample. told (numpy.ndarray) – The sample times from the killed trajectory. xold (numpy.ndarray) – The killed trajectory.
	**kwargs – Keyword arguments, forwarded to `simul_trajectory()`.
Returns:	tnew, xnew – The resampled trajectory.
Return type:	numpy.ndarray, numpy.ndarray

simul_trajectory(x0, t0, **kwargs)¶

Simulate a trajectory until it reaches either set A (score <= 0) or set B (score >= 1).

Keyword Arguments:
Parameters:	x0 (float) – Initial position. t0 (float) – Initial time.
	dt (float) – The time step.
Returns:	t, x – The simulated trajectory.
Return type:	numpy.ndarray, numpy.ndarray

initialize_ensemble(ntraj, **kwargs)¶

Generate the initial ensemble.

Keyword Arguments:
Parameters:	ntraj (int) – Number of trajectories in the ensemble.
	**kwargs – Keyword arguments forwarded to `simul_trajectory()`.

static selectionstep(levels, npart=1)¶

Selection step of the AMS algorithm.

Parameters:	levels (numpy.ndarray) – The list of levels reached by the ensemble members. npart (int, optional) – The number of levels to select. npart=1 corresponds to the last particle method. Note that one level can correspond to several trajectories in the ensemble.
Returns:	killed, survivors – The indices of the killed and surviving ensemble members.
Return type:	numpy.ndarray, numpy.ndarray

Notes

Return the trajectories in the ensemble which performed worse, i.e. the trajectories for which the maximum value of the score function over the trajectory is minimum.

mutationstep(killed_pool, survivor_pool, **kwargs)¶

Mutation step for the AMS algorithm.

Keyword Arguments:
Parameters:	killed_pool (array_like) – The indices of the ensemble members to kill. survivor_pool (array_like) – The indices of the ensemble members to keep.
	**kwargs – Keyword arguments forwarded to `resample()`.

Notes

This is the only method which modifies the state of the ensemble (the AMS object).

run_iter(ntraj, niter, **kwargs)¶

Generate trajectories with the AMS algorithm.

Parameters:	ntraj (int) – The number of trajectories in the initial ensemble. niter (int) – The number of iterations of the algorithm. Arguments (Keywords) – ------------------ – **kwargs – Keyword arguments passed to the “trajectory” method of the dynamics object.
Yields:	trajectory, weight (numpy.ndarray, float) – The generator yields (trajectory, weight) pairs which allows to compute easily the probability associated to each sampled trajectory.

Notes

This method yields first the killed trajectories as the algorithm is iterated, then the trajectories in the final ensemble.

run_resamp(ntraj, niter, **kwargs)¶

Generate trajectories with the AMS algorithm.

Parameters:	ntraj (int) – The number of trajectories in the initial ensemble. niter (int) – The number of iterations of the algorithm. Arguments (Keywords) – ------------------ – **kwargs – Keyword arguments passed to the “trajectory” method of the dynamics object.
Yields:	trajectory, weight (numpy.ndarray, float) – The generator yields (trajectory, weight) pairs which allows to compute easily the probability associated to each sampled trajectory.

Notes

This method yields first the trajectories in the initial ensemble, then the resampled trajectories as the algorithm is iterated.

run_level(ntraj, target_lev, **kwargs)¶

Generate trajectories with the AMS algorithm.

Parameters:	ntraj (int) – The number of trajectories in the initial ensemble. target_lev (float) – The target level. Arguments (Keywords) – ------------------ – **kwargs – Keyword arguments passed to the “trajectory” method of the dynamics object.
Yields:	trajectory, weight (numpy.ndarray, float) – The generator yields (trajectory, weight) pairs which allows to compute easily the probability associated to each sampled trajectory.

Notes

This method yields first the killed trajectories as the algorithm is iterated, then the trajectories in the final ensemble.

class stochrare.rare.ams.TAMS(model, scorefun, duration, **kwargs)¶

Bases: stochrare.rare.ams.AMS

Implement the TAMS algorithm [5].

Parameters:	dynamics (stochrare.dynamics.StochModel object (or a subclass of it)) – The dynamical model; so far we are restricted to SDEs of the form \(dX_t = F(X_t, t) + \sqrt{2D}dW_t\) We only use the trajectory method of the dynamics object. score (function with two arguments.) – The score function \(\xi(t, x)\). duration (float) – The fixed duration for each trajectory.

Notes

This implementation keeps all the information in memory: this should not be suitable for complex dynamics. Similarly, the algorithm is not parallelized, even if the dynamics itself may be.

References

[5]	(1, 2) Lestang, F. Ragone, C.-E. Brehier, C. Herbert and F. Bouchet, J. Stat. Mech. (2018)

resample(time, pos, told, xold, **kwargs)¶

Resample a killed trajectory after a given time.

Keyword Arguments:
Parameters:	time (float) – The time from which to resample. pos (float) – The position from which to resample. told (numpy.ndarray) – The sample times from the killed trajectory. xold (numpy.ndarray) – The killed trajectory.
	**kwargs – Keyword arguments, forwarded to `simul_trajectory()`.
Returns:	tnew, xnew – The resampled trajectory.
Return type:	numpy.ndarray, numpy.ndarray

simul_trajectory(x0, t0, **kwargs)¶

Simulate a trajectory with given initial conditions for a fixed duration.

Keyword Arguments:
Parameters:	x0 (float) – Initial position. t0 (float) – Initial time.
	T (float) – The duration. dt (float) – The time step.
Returns:	t, x – The simulated trajectory.
Return type:	numpy.ndarray, numpy.ndarray

average(ntraj, niter, observable, **kwargs)¶

Estimate the average of an observable using AMS sampling.

Keyword Arguments:
Parameters:	ntraj (int) – The number of initial trajectories in the ensemble. niter (int) – The number of iterations of the AMS algorithm. observable (function with two arguments) – A function of the form O(t, x), where t and x are numpy arrays. It should itself return a numpy array. For instance, it could be a time-independent function of the type `lambda t, x: x2` or a functional of the trajectory such as `lambda t, x: np.array([np.max(x2)])` Note that in the latter case it is crucial to convert the scalar to an array.
	method (function) – The method used to sample the trajectories with the AMS algorithm. It can be one of `run_iter()` (default) or `run_resamp()`. condition (function) – A predicate for conditional averaging. It should be of the form `pred((t,x))` in (True, False).
Returns:	obs – The expectation value of the observable.
Return type:	numpy.ndarray

returntimes(ntraj, niter, **kwargs)¶

Estimate the return time of an observable using AMS sampling.

Keyword Arguments:
Parameters:	ntraj (int) – The number of initial trajectories in the ensemble. niter (int) – The number of iterations of the AMS algorithm.
	method (function) – The method used to sample the trajectories with the AMS algorithm. It can be one of `run_iter()` (default) or `run_resamp()`. observable (function) – The time-dependent observable O(t, x). The default is the score function.
Returns:	a, r(a) – The amplitude and associated return time using the generalized block-maximum method.
Return type:	numpy.ndarray, numpy.ndarray

Classes

`AMS`(model, scorefun[, initcond])	Original version of the Adaptive Multilevel Splitting Algorithm.
`TAMS`(model, scorefun, duration, **kwargs)	Implement the TAMS algorithm [5].