Double Well Spatial Decorrelation
---------------------------------
 
.. code:: ipython2

    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')


.. parsed-literal::

    Populating the interactive namespace from numpy and matplotlib


.. parsed-literal::

    /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"


.. code:: ipython2

    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

.. code:: ipython2

    P = np.array([[0.99, 0.01], 
                  [0.01, 0.99]]);
    T = 50000
    means = [np.array([-1,1]), np.array([1,-1])];
    widths = [np.array([0.3,2]),np.array([0.3,2])];

.. code:: ipython2

    # 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]    

.. code:: ipython2

    dtraj.shape




.. parsed-literal::

    (50000,)



.. code:: ipython2

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



.. image:: doublewellsd_files/doublewellsd_5_0.png


.. code:: ipython2

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




.. parsed-literal::

    <matplotlib.collections.PathCollection at 0x10ca88110>




.. image:: doublewellsd_files/doublewellsd_6_1.png


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.   

.. code:: ipython2

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


.. parsed-literal::

    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

.. code:: ipython2

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


.. parsed-literal::

    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)


.. code:: ipython2

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

.. code:: ipython2

    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*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')
    





.. parsed-literal::

    <matplotlib.text.Text at 0x112a87c50>




.. image:: doublewellsd_files/doublewellsd_12_1.png


.. code:: ipython2

    YSD4 = W.dot(Y); 

.. code:: ipython2

    hist(3*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()




.. parsed-literal::

    <matplotlib.legend.Legend at 0x10cc15b50>




.. image:: doublewellsd_files/doublewellsd_14_1.png