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 [ ]: