Testing the non-Markovian Path Analysis Package



In [1]:

    
import sys
sys.path.append("../nmpath/")
from tools_for_notebook0 import *
%matplotlib inline
from mappers import rectilinear_mapper
from ensembles import Ensemble, DiscreteEnsemble, PathEnsemble, DiscretePathEnsemble

2D Toy model



In [2]:

    
plot_traj([],[])

MC simulation



In [3]:

    
#Generating MC trajectories
mc_traj1_2d = mc_simulation2D(100000)
mc_traj2_2d = mc_simulation2D(10000)

1 - Ensemble class (analysis of continuos trajectories)

Stores an esemble (list) of trajectories (np.arrays). The ensemble could have any number of trajectories including no trajectories at all.

Creating an Ensemble



In [4]:

    
# Empty ensemble with no trajectories
my_ensemble = Ensemble()

From a single trajectory:



In [5]:

    
# from a single trajectory
my_ensemble = Ensemble([mc_traj1_2d],verbose=True)









    



Read 1 (2-dimensional) trajectories of average length 100000.0.

From a list of trajectories:



In [6]:

    
# We have to set list_of_trajs = True

my_list_of_trajs = [mc_traj1_2d, mc_traj2_2d]

my_ensemble = Ensemble(my_list_of_trajs, verbose=True)









    



Read 2 (2-dimensional) trajectories of average length 55000.0.

Ensembles are iterable objects



In [7]:

    
for traj in my_ensemble:
    print(len(traj))

Adding trajectories to the Ensemble

New trajectories can be added to the ensemble as long as there is consistency in the number of variables.



In [8]:

    
my_ensemble = Ensemble(verbose=True)

my_ensemble.add_trajectory(mc_traj1_2d)

my_ensemble.add_trajectory(mc_traj2_2d)









    



Empty ensemble generated

Continuous, Ensemble with 1 (2-dimensional) trajectories
Total number of snapshots: 100000

Continuous, Ensemble with 2 (2-dimensional) trajectories
Total number of snapshots: 110000

"Printing" the ensemble



In [9]:

    
print(my_ensemble)









    



Continuous, Ensemble with 2 (2-dimensional) trajectories
Total number of snapshots: 110000

Defining states and computing MFPTs

The states are considered intervals in the is the class is Ensemble



In [10]:

    
stateA = [[0,pi],[0,pi]]
stateB = [[5*pi,6*pi],[5*pi,6*pi]]

my_ensemble.empirical_mfpts(stateA, stateB)









    



Number of A->B/B->A  events: 32/31






    Out[10]:





{'mfptAB': 1446.125,
 'mfptBA': 2025.7096774193549,
 'std_err_mfptAB': 228.24434699153832,
 'std_err_mfptBA': 368.37653561932768}

Sum of ensembles (ensemble + ensemble)



In [11]:

    
seq1 = mc_simulation2D(20000)
seq2 = mc_simulation2D(20000)

my_e1 = Ensemble([seq1])
my_e2 = Ensemble([seq2])

ensemble1 = my_e1 + my_e2

Another simple example



In [12]:

    
e1 = Ensemble([[1.,2.,3.,4.]],verbose=True)
e2 = Ensemble([[2,3,4,5]])
e3 = Ensemble([[2,1,1,4]])

my_ensembles = [e1, e2, e3]

ensemble_tot = Ensemble([])

for ens in my_ensembles:
    ensemble_tot += ens

#ensemble_tot.mfpts([1,1],[4,4])









    



Read 1 (1-dimensional) trajectories of average length 4.0.

Computing the count matrix and transition matrix



In [13]:

    
n_states = N**2
bin_bounds = [[i*pi for i in range(7)],[i*pi for i in range(7)]]
C1 = my_ensemble._count_matrix(n_states, map_function=rectilinear_mapper(bin_bounds))
print(C1)









    



[[ 4623.    92.     0. ...,     0.     0.     0.]
 [   89.   637.    68. ...,     0.     0.     0.]
 [    0.    75.  3870. ...,     0.     0.     0.]
 ..., 
 [    0.     0.     0. ...,  3594.    83.     0.]
 [    0.     0.     0. ...,    81.   622.    85.]
 [    0.     0.     0. ...,     0.    82.  5410.]]



In [14]:

    
K1 = my_ensemble._mle_transition_matrix(n_states, map_function=rectilinear_mapper(bin_bounds))
print(K1)









    



[[ 0.96072319  0.01911887  0.         ...,  0.          0.          0.        ]
 [ 0.10520095  0.75295508  0.08037825 ...,  0.          0.          0.        ]
 [ 0.          0.0169837   0.8763587  ...,  0.          0.          0.        ]
 ..., 
 [ 0.          0.          0.         ...,  0.86560694  0.01999037  0.        ]
 [ 0.          0.          0.         ...,  0.0971223   0.74580336
   0.10191847]
 [ 0.          0.          0.         ...,  0.          0.01459854
   0.96314759]]

2 - PathEnsemble class

Creating a path ensemble object



In [15]:

    
#p_ensemble = PathEnsemble()

From ensemble



In [16]:

    
p_ensemble = PathEnsemble.from_ensemble(my_ensemble, stateA, stateB)

print(p_ensemble)









    



Continuous, PathEnsemble with 32 (2-dimensional) trajectories
Total number of snapshots: 46308

MFPTs



In [17]:

    
p_ensemble.empirical_mfpts(stateA, stateB)









    



WARNING: No B->A events observed
Number of A->B/B->A  events: 32/0






    Out[17]:





{'mfptAB': 1446.125,
 'mfptBA': 'NaN',
 'std_err_mfptAB': 228.24434699153832,
 'std_err_mfptBA': 'NaN'}

Count matrix



In [18]:

    
print(p_ensemble._count_matrix(n_states, mapping_function2D))









    



[[ 4449.    90.     0. ...,     0.     0.     0.]
 [   73.   488.    48. ...,     0.     0.     0.]
 [    0.    41.  1738. ...,     0.     0.     0.]
 ..., 
 [    0.     0.     0. ...,  1526.    27.     0.]
 [    0.     0.     0. ...,    18.   129.    14.]
 [    0.     0.     0. ...,     0.     0.     0.]]



In [19]:

    
#clusters = p_ensemble.cluster(distance_metric = 'RMSD', n_cluster=10, method = 'K-means')

3 - DiscreteEnsemble class

We can generate a discrete ensemble from the same mapping function and we should obtain exaclty the same result:



In [20]:

    
d_ens = DiscreteEnsemble.from_ensemble(my_ensemble, mapping_function2D)
print(d_ens)









    



Discrete, DiscreteEnsemble with 2 (1-dimensional) trajectories
Total number of snapshots: 110000

Count matrix and transition matrix



In [21]:

    
C2 = d_ens._count_matrix(n_states)
print(C2)









    



[[ 4623.    92.     0. ...,     0.     0.     0.]
 [   89.   637.    68. ...,     0.     0.     0.]
 [    0.    75.  3870. ...,     0.     0.     0.]
 ..., 
 [    0.     0.     0. ...,  3594.    83.     0.]
 [    0.     0.     0. ...,    81.   622.    85.]
 [    0.     0.     0. ...,     0.    82.  5410.]]



In [22]:

    
K2= d_ens._mle_transition_matrix(n_states)
print(K2)









    



[[ 0.96072319  0.01911887  0.         ...,  0.          0.          0.        ]
 [ 0.10520095  0.75295508  0.08037825 ...,  0.          0.          0.        ]
 [ 0.          0.0169837   0.8763587  ...,  0.          0.          0.        ]
 ..., 
 [ 0.          0.          0.         ...,  0.86560694  0.01999037  0.        ]
 [ 0.          0.          0.         ...,  0.0971223   0.74580336
   0.10191847]
 [ 0.          0.          0.         ...,  0.          0.01459854
   0.96314759]]

Defining states and computing MFPTs

The states are now considered sets, defining the states as follow we should obtain the same results



In [23]:

    
stateA = [0]
stateB = [N*N-1]

d_ens.empirical_mfpts(stateA, stateB)









    



Number of A->B/B->A  events: 32/31






    Out[23]:





{'mfptAB': 1446.125,
 'mfptBA': 2025.7096774193549,
 'std_err_mfptAB': 228.24434699153832,
 'std_err_mfptBA': 368.37653561932768}

Generating a Discrete Ensemble from the transition matrix



In [24]:

    
d_ens2 = DiscreteEnsemble.from_transition_matrix(K2, sim_length = 100000)



In [25]:

    
#d_ens2.mfpts(stateA,stateB)

4 - DiscretePathEnsemble class

Creating the DPE

From Ensemble



In [26]:

    
dpathEnsemble = DiscretePathEnsemble.from_ensemble(my_ensemble, stateA, stateB, mapping_function2D)
print(dpathEnsemble)









    



Discrete, DiscretePathEnsemble with 32 (1-dimensional) trajectories
Total number of snapshots: 46308



In [27]:

    
#MFPT from the transition matrix
dpathEnsemble.nm_mfpt(ini_probs = None, n_states = N*N)









    Out[27]:





1446.1249999999975

From the transition matrix



In [28]:

    
n_paths = 200

dpathEnsemble2 = DiscretePathEnsemble.from_transition_matrix\
                (K2, stateA = stateA, stateB = stateB, n_paths = n_paths,ini_pops = [1])
    
print(dpathEnsemble2)









    



Discrete, DiscretePathEnsemble with 200 (1-dimensional) trajectories
Total number of snapshots: 198959

Fundamental sequence



In [29]:

    
FSs, weights, count = dpathEnsemble2.weighted_fundamental_sequences(K2)
size = len(FSs)

paths = dpathEnsemble2.trajectories
print(count)

Plotting paths A -> B



In [30]:

    
discrete = [True for i in range(size)]

plot_traj([[paths[i],[]] for i in range(size)] , discrete, \
          line_width=0.2, std=0.5, color='k', title = '{} paths A->B'.format(n_paths))

Plotting Fundamental Sequences A -> B



In [31]:

    
plot_traj([[FSs[i],[]] for i in range(size)] ,discrete, \
          line_width=0.5, std=0.2, color='k', title = '{} FSs A->B'.format(n_paths))



In [32]:

    
lw = [weights[i]*100 for i in range(size)]



In [33]:

    
#np.random.seed(12)
plot_traj([[FSs[i],[]] for i in range(size)] ,discrete=[True for i in range(size)],\
          line_width = lw,std = 0.002, alpha=0.25)



In [ ]: