4.3.1. Markov Chain Monte Carlo on hopping graph — hop.MCMC

The hop.MCMC module uses the information encoded in a hopping graph to set up a Markov Chain Monte Carlo sampling procedure that allows one to rapidly generate site occupancy distributions that are distributed in the same way as the one sampled from MD.

The most convenient entry point is the hop.MCMC.run() function

M = MCMC.run(filename='hopgraph.pickle',Ntotal=<int>)

It takes as input a stored hopgraph and immediately runs an MCMC run of Ntotal steps. The output is a MCMCsampler object. It contains the ‘trajectory’ and useful analysis functions. (Use interactive introspection in ipython to explore the possibilities of the object.)

Notes on the algorithm:

• some sort of dynamic lattice Monte Carlo with very simple acceptance probabilities (0 or 1, if there’s no space on the site, and 1 if there is)

  … is ‘MCMC’ actually the proper description?

• Extension to multiply occupied sites: use the site occupancy distributions from siteanalysis, and replace the unconditional move by an acceptance probability == s_i(n)

• I am currently using time-forward (out-of) and time-backward (into) site moves (the latter inspired by coupling from the past).

4.3.1.1. Classes and functions

class hop.MCMC.MCMCsampler(h=None, min_hops_observed=1, filename=None)

Generate an equilibrium distribution of states from a hop graph.

Initialize the Markov Chain Monte Carlo sample with a HoppingGraph.

M = MCMCsampler(HoppingGraph)

Build a Markov Chain model from a <HoppingGraph> (with edges deleted that have less than <min_hops_observed> hops).

autocorrelation(start=None, stop=None, step=None, **kwargs)

Calculates the auto correlation function for all site trajectories.

averaged_autocorrelation(step=None, **kwargs)

Calculates the ACF or each site by resampling from the whole trajectory.

mean(acf), standardev(acf) = averaged_autocorrelation(**kwargs)

Arguments:

step only take every <step> from the trajectory (None == 1)

??? step > 1 seems to take LONGER ???

length length (in frames) of the ACF (default: 1/2*len(series)) sliding_window repeat ACF calculation every N frames (default: len(series)/100)

Returns:

mean_acf average over all resampled acfs per site, shape = (Nsites,length) std_acf standard deviation or the resampled acfs, shape = (Nsites,length)

See also for kwargs:

firstsiteindex

State array index of the first site after bulk.

index2node

Translates sequential array index to node label (in graph).

init_state(Nbulk=10000.0)

Initialize state with 1 particle per standard site and Nbulk for the bulk site.

mean()

Mean for each site (excl bulk).

mean_std()

Returns site labels, mean, and standard deviation for each site (excl bulk).

node2index

Translates node label (in graph) to the sequential array index.

occupancy()

Ensemble averaged occupancy (over ALL states incl bulk) and fluctuation.

occupancy_mean_correl()

Calculates the correlation coefficient between simulation and MCMC occupancies.

occupancy_std_correl()

Calculates the correlation coefficient between simulation and MCMC occupancy fluctuations.

plot(filename=None, plot_skip=None)

Plot density plot of the saved configurations in states[].

plot_correl(legend=True, **kwargs)

Plot the occupancy from the MD simulation vs the MCMC one.

plot_occupancy(legend=True, **kwargs)

Plot site label vs <N> +/- std(N).

legend True: add legend, False: return (line,description) **kwargs additional arguments for errbar plot such as color=’k’, fmt=’o’

run(Ntotal=500000, Nskip=1000, verbosity=None)

MCMC run multiple cycles of lebgth <Nskip> scans for a total of <Ntotal>.

run(Ntotal=500000,Nskip=1000)

Starts from the current configuration in state. Creates the collection of configurations states: one state every Nskip steps

sample(max_iter=10000, record_iterations=True)

Run <max_iter> Monte Carlo site moves.

sample(max_iter=10000)

Runs a batch of MCMC moves.

sites

Translates sequential array index to node label (in graph).

statevector

State as a numpy array; the corresponding nodes are state.keys()

std()

Standard deviation for each site.

class hop.MCMC.MultiPscan(repeat=10, **pscanargs)

Run Pscan(**pscanargs) <repeat> times and collect all Pscan objects in list.

pscans = MultiPscan(repeat=10, parameter=’Ntotal’, pvalues=[1e4,2.5e4,….], …) See Pscan() for description of pscanargs.

class hop.MCMC.Pscan(parameter, pvalues=None, filename='hopgraph.pickle', Ntotal=1000000.0, **kwargs)

Run a MCMC sampler for a number of parameter values.

Sample on hopping graph for different values of <parameter> = p.

P = Pscan(parameter=<string>,pvalues=<sequence>,filename=<filename>,**kwargs)

<parameter> must be a keyword argument to hop.MCMC.run(); the parameter overrides any default values that may have been set. For instance, <parameter> can be ‘Ntotal’ or ‘filename’.

kwargs: all other kwargs are directly passed on to MCMC.run().

occupancy_mean_correl()

Returns X=pvalues, Y=occupancy_mean_correlations.

plot_occupancy(**kwargs)

Plot <n_i> (site occupancy from MCMC) for all parameter values.

(See _plotter())

plot_occupancy_mean_correl(**kwargs)

Plot MD occupancy vs MCMC occupancy.

plot_correl(colorscale=’linear’|’log’)

(See _plotter())

plot_states(maxcolumns=2)

Plot all state ‘trajectories’ as a tiled plot.

save(filename='pscan.pickle')

Save pscan object to pickle file.

save(pscan.pickle)

Load with

import cPickle myPscan = cPickle.load(open(‘pscan.pickle’))

hop.MCMC.multi_plot(plist, plottype='whisker', Nequil=10000, funcname='occupancy_mean_correl', **kwargs)

Display a collection of functions.

multi_plot(plist,plottype=’whisker’,Nequil=10000,funcname=’occupancy_mean_correl’,**kwargs)

The function is obtained from a method call on the objects in plist. The assumption is that these are functions of Ntotal (if not, set Nequil=0; Nequil is added to x). Each object is a different realization, e.g. multiple MCMC runs.

plottype ‘whisker’ (whisker plot), ‘standard’ (average and standard deviations) Nequil correction, added to x funcname string; a method of the objects in plist that does EXACTLY the following:

x,y = obj.funcname() where x and y are numpy arrays of equal length

**kwargs color, boxcolor, mediancolor, capsize

hop.MCMC.run(filename='hopgraph.pickle', Ntotal=500000, Nskip=1000, Nequil=10000)

Perform Markov Chain Monte Carlo on a model derived from the hopping graph.