Triple Well Spatial Decorrelation

import numpy as np
import pyemma.msm as msm
import msmtools.generation as msmgen
import msmtools.analysis as msmana
import pyemma.coordinates as coor
import matplotlib.pylab as plt
import anca
%pylab inline
plt.style.use('ggplot')

Populating the interactive namespace from numpy and matplotlib

/Users/fxp/anaconda2/lib/python2.7/site-packages/IPython/core/magics/pylab.py:161: UserWarning: pylab import has clobbered these variables: ['plt']
%matplotlib prevents importing * from pylab and numpy
  "n`%matplotlib` prevents importing * from pylab and numpy"

def assign(X, cc):
    T = X.shape[0]
    I = np.zeros((T),dtype=int)
    for t in range(T):
        dists = X[t] - cc
        dists = dists ** 2
        I[t] = np.argmin(dists)
    return I

P = np.array([[0.75, 0.05, 0.2],
              [0.05, 0.75, 0.2],
              [0.2, 0.05, 0.75]]);
T = 100000
means = [np.array([-1,1]), np.array([1,-1]), np.array([2.5,2.5])];
widths = [np.array([0.3,2]),np.array([0.3,2]), np.array([0.3,2])];

# continuous trajectory
X = np.zeros((T, 2))
# hidden trajectory
dtraj = msmgen.generate_traj(P, T)
for t in range(T):
    s = dtraj[t]
    X[t,0] = widths[s][0] * numpy.random.randn() + means[s][0]
    X[t,1] = widths[s][1] * numpy.random.randn() + means[s][1]

dtraj.shape

(100000,)

plt.plot(dtraj[0:500]);

plt.figure(figsize=(4,7))
plt.scatter(X[:,0], X[:,1], marker = 'o', color=[0.6,0.6,0.6])

<matplotlib.collections.PathCollection at 0x10c93e550>

Spatial Decorrelation of Order 2 (SD2)

Parameters:

data – a 3n x T data matrix (number 3 is due to the x,y,z coordinates for each atom). Maybe a numpy
array or a matrix where,

n: size of the protein

T: number of snapshots of MD trajectory

m – dimensionality of the subspace we are interested in; Default value is None, in which case m = n
verbose – print information on progress. Default is true.

Returns:

A 3n x m matrix U (NumPy matrix type), such that Y = U * data is a 2nd order spatially whitened coordinates extracted from the 3n x T data matrix. If m is omitted, U                is a square 3n x 3n matrix.

from pyANCA import SD2
(Y, S, B, U) = SD2.SD2(X, m=2);

2nd order Spatial Decorrelation -> Looking for 2 sources
2nd order Spatial Decorrelation -> Removing the mean value
2nd order Spatial Decorrelation -> Whitening the data

Spatial Decorrelation Module of Order 4 (SD4)

Parameters:

Y -- an mxT spatially whitened matrix (m dimensionality of subspace, T snapshots).
May be a numpy array or a matrix where

m -- dimensionality of the subspace we are interested in. Defaults to None, in which case m=n.
T -- number of snapshots of MD trajectory

U -- whitening matrix obtained after doing the PCA analysis on m components of real data

verbose -- print info on progress. Default is True.

Returns:

W -- a separating matrix for spatial decorrelation of order 4

from pyANCA import SD4
W = SD4.SD4(Y, m=2, U=U)

4th order Spatial Decorrelation -> Estimating cumulant matrices
SD4 -> Contrast optimization by joint diagonalization
SD4 -> Sweep #  0 completed in 1 rotations
SD4 -> Sweep #  1 completed in 0 rotations
SD4 -> Total of 1 Givens rotations
SD4 -> Sorting the components
SD4 -> Fixing the signs
(2, 2)

def draw_arrow(a, v, color):
    plt.arrow(0, 0, a*v[0], a*v[1], color=color, width=0.02, linewidth=3)

plt.figure(figsize=(4,7))
scatter(X[:,0], X[:,1], marker = 'o', color=[0.6,0.6,0.6])

plt.arrow(0, 0, 7*U[0,0], 12*U[0,1], color='red', width=0.02, linewidth=3);
plt.text(-0.0, 6.5, 'SD2', color='red', fontsize=20, fontweight='bold', rotation='horizontal')

#plt.arrow(0, 0, 3*V[0,0], 4*V[0,1], color='blue', width=0.02, linewidth=3);
#plt.text(-1.5, 3.5, 'TD2', color='blue', fontsize = 20, fontweight='bold', rotation='horizontal')

plt.arrow(0, 0, 7*W[0,0], 9*W[0,1], color='orange', width=0.02, linewidth=3);
plt.text(1.5, 3.5, 'SD4', color='orange', fontsize=20, fontweight='bold', rotation='horizontal')

<matplotlib.text.Text at 0x1126c3a90>

YSD4 = W.dot(Y);

hist(Y[0,:].T, bins=50, histtype='step', linewidth=3, label='SD2', color='blue')
hist(4*YSD4[0,:].T, bins=50, histtype='step', linewidth=3, label='SD4', color='red')

xlabel('essential coordinate (1st principal or independent component)')
ylabel('projected histogram')
legend()

<matplotlib.legend.Legend at 0x111e0f7d0>