# Single Replica TIS on a Toy Model

This continues the work from the toy MSTIS example. Taking the results from that, we can create a fixed bias single replica TIS simulation.

Normally you would iterate the fixed bias calculation several times in order to converge toward a flat histogram. In this particular example, we start with a well-converged result from the previous example, so we don’t need to iterate.

toy_mstis_5_srtis

# Single Replica TIS¶

This notebook shows how to run single replica TIS move scheme. This assumes you can load engine, network, and initial sample from a previous calculation.

In [ ]:
%matplotlib inline
import openpathsampling as paths
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

In [ ]:
from openpathsampling.visualize import PathTreeBuilder, PathTreeBuilder
from IPython.display import SVG, HTML

def ipynb_visualize(movevis):
"""Default settings to show a movevis in an ipynb."""
view = movevis.renderer
view.zoom = 1.5
view.scale_y = 18
view.scale_th = 20
view.font_size = 0.4
return view


Open the storage and load things from it.

In [ ]:
old_store = paths.AnalysisStorage("mstis.nc")
#old_store = paths.Storage("mstis.nc")  # if not actually doing analysis, but loading network, etc

In [ ]:
network = old_store.networks[0]
engine = old_store.engines[0]
template = old_store.snapshots[0]


One of the points of SRTIS is that we use a bias (which comes from an estimate of the crossing probability) in order to improve our sampling.

In [ ]:
# this is how we would get it out of a simulation (although the actual simulation here has bad stats)
# first, we need the crossing probabilities, which we get when we calculate the rate
network.hist_args['max_lambda'] = { 'bin_width' : 0.02, 'bin_range' : (0.0, 0.5) }
network.hist_args['pathlength'] = { 'bin_width' : 5, 'bin_range' : (0, 150) }
rates = network.rate_matrix(old_store.steps)

In [ ]:
# just use the analyzed network to make the bias
bias = paths.SRTISBiasFromNetwork(network)
bias.df

In [ ]:
# For better stats, use the results that I got from a 20k MC step run
# We can create fake TCPs and force them on the network.

tcp_A = paths.numerics.LookupFunction.from_dict(
{0.2: 1.0,
0.3: 0.13293104100673198,
0.4: 0.044370838092911397,
0.5: 0.021975696374764188}
)
tcp_B = paths.numerics.LookupFunction.from_dict(
{0.2: 1.0,
0.3: 0.13293104100673198,
0.4: 0.044370838092911397,
0.5: 0.021975696374764188}
)
tcp_C = paths.numerics.LookupFunction.from_dict(
{0.2: 1.0,
0.3: 0.19485705066078274,
0.4: 0.053373003923696649,
0.5: 0.029175949467020165}
)

# load states for identification purposes
stateA = old_store.volumes['A']
stateB = old_store.volumes['B']
stateC = old_store.volumes['C']

# use the sampling transitions; in MSTIS, these are also stored in from_state
network.from_state[stateA].tcp = tcp_A
network.from_state[stateB].tcp = tcp_B
network.from_state[stateC].tcp = tcp_C

In [ ]:
bias = paths.SRTISBiasFromNetwork(network)
bias.df


Here we actually set up the SRTIS move scheme for the given network. It only requires one line:

In [ ]:
scheme = paths.SRTISScheme(network, bias=bias, engine=engine)


Now we'll visualize the SRTIS move scheme.

In [ ]:
movevis = paths.visualize.MoveTreeBuilder()
#movevis.mover(scheme.move_decision_tree(), network.all_ensembles)
#SVG(ipynb_visualize(movevis).to_svg())


Next we need to set up an appropriate single-replica initial sampleset. We'll take the last version of from one of the outer TIS ensembles.

In [ ]:
final_samp0 = old_store.steps[len(old_store.steps)-1].active[network.sampling_ensembles[-1]]

In [ ]:
sset = paths.SampleSet([final_samp0])


Finally, we set up the new storage file and the new simulation.

In [ ]:
storage = paths.Storage("srtis.nc", "w")
storage.save(template)

In [ ]:
srtis = paths.PathSampling(
storage=storage,
sample_set=sset,
move_scheme=scheme
)

In [ ]:
n_steps_to_run = int(scheme.n_steps_for_trials(
mover=scheme.movers['minus'][0],
n_attempts=1
))
print(n_steps_to_run)

In [ ]:
# logging creates ops_output.log file with details of what the calculation is doing
#import logging.config
#logging.config.fileConfig("logging.conf", disable_existing_loggers=False)

In [ ]:
%%time
multiplier = 2
srtis.run_until(multiplier*n_steps_to_run)

In [ ]:
#storage.close()


From here, we'll be doing the analysis of the SRTIS run.

In [ ]:
%%time
#storage = paths.AnalysisStorage("srtis.nc")
#scheme = storage.schemes[0]

In [ ]:
scheme.move_summary(storage.steps)

In [ ]:
scheme.move_summary(storage.steps, 'shooting')

In [ ]:
scheme.move_summary(storage.steps, 'minus')

In [ ]:
scheme.move_summary(storage.steps, 'repex')

In [ ]:
scheme.move_summary(storage.steps, 'pathreversal')

In [ ]:
replica = storage.samplesets[0].samples[0].replica
ensemble_trace = paths.trace_ensembles_for_replica(replica, storage.steps)
print len(ensemble_trace)

In [ ]:
srtis_ensembles = scheme.network.sampling_ensembles+scheme.network.special_ensembles['ms_outer'].keys()
srtis_ensemble_numbers = {e : srtis_ensembles.index(e) for e in srtis_ensembles}
# this next is just for pretty printing
srtis_numbers_ensemble = {srtis_ensemble_numbers[e] : e for e in srtis_ensemble_numbers}
for k in sorted(srtis_numbers_ensemble.keys()):
print k, ":", srtis_numbers_ensemble[k].name

In [ ]:
plt.plot([srtis_ensemble_numbers[e] for e in ensemble_trace])

In [ ]:
count = 0
for i in range(len(ensemble_trace)-1):
[this_val, next_val] = [srtis_ensemble_numbers[ensemble_trace[k]] for k in [i,i+1]]
if this_val == 1 and next_val == 0:
count += 1
count

In [ ]:
hist_numbers = [srtis_ensemble_numbers[e] for e in ensemble_trace]
bins = [i-0.5 for i in range(len(srtis_ensembles)+1)]

In [ ]:
plt.hist(hist_numbers, bins=bins);

In [ ]:
import pandas as pd
hist = paths.analysis.Histogram(bin_width=1.0, bin_range=[-0.5,9.5])
colnames = {i : srtis_numbers_ensemble[i].name for i in range(len(srtis_ensembles))}
df = pd.DataFrame(columns=[colnames[i] for i in colnames])

In [ ]:
for i in range(len(hist_numbers)):

plt.pcolormesh(df.fillna(0.0), cmap="bwr", vmin=0.0, vmax=0.2);