In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import openpathsampling as paths
import openpathsampling.storage as st

storage = st.AnalysisStorage("mstis.nc")

In [2]:
len(storage.steps)


Out[2]:
501

In [3]:
mstis = storage.networks.load(0)
retis = mstis.sampling_transitions[0]
scheme = storage.schemes[0]

Defining labels

Everything is actually done in terms of ensembles. We can map the ensembles to any labels. In our case, we use the initial replica ID associated with the ensemble. We use this as both the numeric and string label.


In [4]:
sset0 = storage.samplesets[0]
numeric_labels = { s.ensemble : s.replica for s in sset0}
string_labels = { s.ensemble : str(s.replica) for s in sset0 }
numeric_to_string = { numeric_labels[e] : string_labels[e] for e in numeric_labels.keys()}

Trace of ensemble visited by a replica

In the plot below, you'll see we set the labels on the axis as sset0[e].replica, so we use the replica ID associated with the ensemble in the first timestep.


In [5]:
%%time
trace_1 = paths.trace_ensembles_for_replica(0, storage.steps)


CPU times: user 6.79 ms, sys: 2.02 ms, total: 8.81 ms
Wall time: 7.5 ms

In [6]:
plt.plot([numeric_labels[e] for e in trace_1])


Out[6]:
[<matplotlib.lines.Line2D at 0x11014ce90>]

Replica flow


In [7]:
repx_net = paths.ReplicaNetwork(scheme, storage.steps)
print repx_net


<openpathsampling.analysis.replica_network.ReplicaNetwork object at 0x11062e190>

In [8]:
flow = repx_net.flow(bottom=retis.minus_ensemble, top=retis.ensembles[-1])
print flow


{<openpathsampling.ensemble.TISEnsemble object at 0x10de2ead0>: 0.0, <openpathsampling.ensemble.TISEnsemble object at 0x10df0ae50>: 0.0, <openpathsampling.ensemble.MinusInterfaceEnsemble object at 0x10dee7310>: 0.0, <openpathsampling.ensemble.TISEnsemble object at 0x10defa890>: 0.8551068883610451, <openpathsampling.ensemble.TISEnsemble object at 0x10dff62d0>: 0.0, <openpathsampling.ensemble.TISEnsemble object at 0x10e022cd0>: 0.0, <openpathsampling.ensemble.MinusInterfaceEnsemble object at 0x10def2610>: 0.0, <openpathsampling.ensemble.TISEnsemble object at 0x10df19110>: 0.0, <openpathsampling.ensemble.UnionEnsemble object at 0x10df14150>: 0.0, <openpathsampling.ensemble.TISEnsemble object at 0x10de2ed90>: 0.0, <openpathsampling.ensemble.TISEnsemble object at 0x10de11990>: 0.0, <openpathsampling.ensemble.MinusInterfaceEnsemble object at 0x10df14bd0>: 1.0, <openpathsampling.ensemble.TISEnsemble object at 0x10de2ee50>: 0.0}

In [9]:
flow_num = {numeric_labels[k] : flow[k] for k in flow.keys()}
print flow_num


{0: 0.0, 1: 0.0, 2: 0.0, 3: 0.0, 4: 0.8551068883610451, 5: 0.0, 6: 0.0, 7: 0.0, 8: 0.0, 9: 0.0, -1: 1.0, -3: 0.0, -2: 0.0}

In [10]:
sorted_vals = []
for k in sorted(flow_num.keys()):
    sorted_vals.append(flow_num[k])

In [11]:
plt.plot(sorted(flow_num.keys()), sorted_vals)


Out[11]:
[<matplotlib.lines.Line2D at 0x110774c50>]

Trips

Now we calculate "up" trips, "down" trips, and round-trips.


In [12]:
repx_net.trips(bottom=retis.minus_ensemble, top=retis.ensembles[-1])


Out[12]:
{'down': [], 'round': [], 'up': []}

Transition matrix

The transition matrix is the (unsymmetrized) matrix of the transition probabilities. By default, it automatically sets its order from the Cuthill-McKee reverse ordering algorithm.


In [13]:
repx_net.transition_matrix()


Out[13]:
12 6 7 8 9 5 4 3 10 11 0 1 2
12 0.000000 0.027397 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
6 0.027397 0.000000 0.082192 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
7 0.000000 0.082192 0.000000 0.047945 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
8 0.000000 0.000000 0.047945 0.000000 0.006849 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
9 0.000000 0.000000 0.000000 0.006849 0.000000 0.020548 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
5 0.000000 0.000000 0.000000 0.000000 0.020548 0.000000 0.054795 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
4 0.000000 0.000000 0.000000 0.000000 0.000000 0.054795 0.000000 0.054795 0.000000 0.000000 0.000000 0.000000 0.000000
3 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.054795 0.000000 0.013699 0.000000 0.000000 0.000000 0.000000
10 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.013699 0.000000 0.000000 0.000000 0.000000 0.000000
11 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.027397 0.000000 0.000000
0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.027397 0.000000 0.041096 0.000000
1 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.041096 0.000000 0.075342
2 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.075342 0.000000

If you would like to set a different order, that can be done by providing a list of the ensembles in whatever order you choose:


In [14]:
import numpy as np
perm = np.random.permutation(len(mstis.all_ensembles))
print perm


[ 5  8 10 11  9  0  4  2 12  7  6  3  1]

In [15]:
order = [mstis.all_ensembles[p] for p in perm]
repx_net.transition_matrix(index_order=order)


Out[15]:
5 8 10 11 9 0 4 2 12 7 6 3 1
5 0.000000 0.000000 0.000000 0.000000 0.020548 0.000000 0.054795 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
8 0.000000 0.000000 0.000000 0.000000 0.006849 0.000000 0.000000 0.000000 0.000000 0.047945 0.000000 0.000000 0.000000
10 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.013699 0.000000
11 0.000000 0.000000 0.000000 0.000000 0.000000 0.027397 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
9 0.020548 0.006849 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
0 0.000000 0.000000 0.000000 0.027397 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.041096
4 0.054795 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.054795 0.000000
2 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.075342
12 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.027397 0.000000 0.000000
7 0.000000 0.047945 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.082192 0.000000 0.000000
6 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.027397 0.082192 0.000000 0.000000 0.000000
3 0.000000 0.000000 0.013699 0.000000 0.000000 0.000000 0.054795 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
1 0.000000 0.000000 0.000000 0.000000 0.000000 0.041096 0.000000 0.075342 0.000000 0.000000 0.000000 0.000000 0.000000

Mixing matrix

Same as the transition matrix $T$, but $\frac{1}{2}(T+T^T)$.


In [16]:
repx_net.mixing_matrix()


Out[16]:
12 6 7 8 9 5 4 3 10 11 0 1 2
12 0.000000 0.027397 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
6 0.027397 0.000000 0.082192 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
7 0.000000 0.082192 0.000000 0.047945 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
8 0.000000 0.000000 0.047945 0.000000 0.006849 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
9 0.000000 0.000000 0.000000 0.006849 0.000000 0.020548 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
5 0.000000 0.000000 0.000000 0.000000 0.020548 0.000000 0.054795 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
4 0.000000 0.000000 0.000000 0.000000 0.000000 0.054795 0.000000 0.054795 0.000000 0.000000 0.000000 0.000000 0.000000
3 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.054795 0.000000 0.013699 0.000000 0.000000 0.000000 0.000000
10 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.013699 0.000000 0.000000 0.000000 0.000000 0.000000
11 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.027397 0.000000 0.000000
0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.027397 0.000000 0.041096 0.000000
1 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.041096 0.000000 0.075342
2 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.075342 0.000000

Making a pretty picture


In [17]:
repxG = paths.ReplicaNetworkGraph(repx_net)

In [18]:
# draw('graphviz') gives better results, but requires pygraphviz
repxG.draw('spring')


Blue is a minus interface, red is a normal interface. Multiple state outer interfaces (not in this example) would be green.

Alternate way of calculating transitions

There's another, perhaps better, way to calculate transitions. This does double count, but doesn't care if the how the transition happened (only that it did).


In [19]:
transitions = repx_net.transitions_from_traces(storage.steps)

In [20]:
for (k1, k2)  in transitions.keys():
    print numeric_labels[k1], numeric_labels[k2], transitions[(k1, k2)]


5 2 8
2 -3 2
9 8 3
-3 2 2
4 0 6
6 7 12
1 0 11
-1 4 4
2 5 8
4 -1 4
8 9 3
0 4 6
9 3 1
3 7 7
-2 6 4
7 6 12
7 3 7
5 8 8
8 5 8
0 1 11
3 9 1
6 -2 4

In [21]:
for (k1, k2) in repx_net.analysis['n_accepted'].keys():
        print numeric_labels[k1], numeric_labels[k2], repx_net.analysis['n_accepted'][(k1, k2)]


5 2 8
2 -3 2
9 8 3
-3 2 2
4 0 6
6 7 12
4 -1 4
-1 4 4
2 5 8
1 0 11
8 9 3
0 4 6
3 9 1
3 7 7
-2 6 4
7 6 12
7 3 7
8 5 8
5 8 8
0 1 11
9 3 1
6 -2 4

In [ ]:


In [ ]: