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.
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)):
hist.add_data_to_histogram([hist_numbers[i]])
if i % 100 == 0:
normalized = hist.normalized()
local_df = pd.DataFrame([normalized.values()], index=[i], columns=[colnames[k] for k in normalized.keys()])
df = df.append(local_df)
In [ ]:
plt.pcolormesh(df.fillna(0.0), cmap="bwr", vmin=0.0, vmax=0.2);
plt.gca().invert_yaxis()
plt.colorbar()
In [ ]: