Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2015-2020 The plumed team
       3             :    (see the PEOPLE file at the root of the distribution for a list of names)
       4             : 
       5             :    See http://www.plumed.org for more information.
       6             : 
       7             :    This file is part of plumed, version 2.
       8             : 
       9             :    plumed is free software: you can redistribute it and/or modify
      10             :    it under the terms of the GNU Lesser General Public License as published by
      11             :    the Free Software Foundation, either version 3 of the License, or
      12             :    (at your option) any later version.
      13             : 
      14             :    plumed is distributed in the hope that it will be useful,
      15             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             :    GNU Lesser General Public License for more details.
      18             : 
      19             :    You should have received a copy of the GNU Lesser General Public License
      20             :    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : #include "core/ActionShortcut.h"
      23             : #include "core/ActionPilot.h"
      24             : #include "core/ActionWithValue.h"
      25             : #include "core/ActionRegister.h"
      26             : #include "core/PlumedMain.h"
      27             : #include "core/ActionSet.h"
      28             : 
      29             : //+PLUMEDOC DIMRED CLASSICAL_MDS
      30             : /*
      31             : Create a low-dimensional projection of a trajectory using the classical multidimensional
      32             :  scaling algorithm.
      33             : 
      34             : Multidimensional scaling (MDS) is similar to what is done when you make a map. You start with distances
      35             : between London, Belfast, Paris and Dublin and then you try to arrange points on a piece of paper so that the (suitably scaled)
      36             : distances between the points in your map representing each of those cities are related to the true distances between the cities.
      37             : Stating this more mathematically MDS endeavors to find an <a href="http://en.wikipedia.org/wiki/Isometry">isometry</a>
      38             : between points distributed in a high-dimensional space and a set of points distributed in a low-dimensional plane.
      39             : In other words, if we have \f$M\f$ \f$D\f$-dimensional points, \f$\mathbf{X}\f$,
      40             : and we can calculate dissimilarities between pairs them, \f$D_{ij}\f$, we can, with an MDS calculation, try to create \f$M\f$ projections,
      41             : \f$\mathbf{x}\f$, of the high dimensionality points in a \f$d\f$-dimensional linear space by trying to arrange the projections so that the
      42             : Euclidean distances between pairs of them, \f$d_{ij}\f$, resemble the dissimilarities between the high dimensional points.  In short we minimize:
      43             : 
      44             : \f[
      45             : \chi^2 = \sum_{i \ne j} \left( D_{ij} - d_{ij} \right)^2
      46             : \f]
      47             : 
      48             : where \f$D_{ij}\f$ is the distance between point \f$X^{i}\f$ and point \f$X^{j}\f$ and \f$d_{ij}\f$ is the distance between the projection
      49             : of \f$X^{i}\f$, \f$x^i\f$, and the projection of \f$X^{j}\f$, \f$x^j\f$.  A tutorial on this approach can be used to analyze simulations
      50             : can be found in the tutorial \ref lugano-5 and in the following <a href="https://www.youtube.com/watch?v=ofC2qz0_9_A&feature=youtu.be" > short video.</a>
      51             : 
      52             : \par Examples
      53             : 
      54             : The following command instructs plumed to construct a classical multidimensional scaling projection of a trajectory.
      55             : The RMSD distance between atoms 1-256 have moved is used to measure the distances in the high-dimensional space.
      56             : 
      57             : \plumedfile
      58             : data: COLLECT_FRAMES ATOMS=1-256
      59             : mat: EUCLIDEAN_DISSIMILARITIES USE_OUTPUT_DATA_FROM=data
      60             : mds: CLASSICAL_MDS USE_OUTPUT_DATA_FROM=mat NLOW_DIM=2
      61             : OUTPUT_ANALYSIS_DATA_TO_COLVAR USE_OUTPUT_DATA_FROM=mds FILE=rmsd-embed
      62             : \endplumedfile
      63             : 
      64             : The following section is for people who are interested in how this method works in detail. A solid understanding of this material is
      65             : not necessary to use MDS.
      66             : 
      67             : \section dim-sec Method of optimization
      68             : 
      69             : The stress function can be minimized using a standard optimization algorithm such as conjugate gradients or steepest descent.
      70             : However, it is more common to do this minimization using a technique known as classical scaling.  Classical scaling works by
      71             : recognizing that each of the distances $D_{ij}$ in the above sum can be written as:
      72             : 
      73             : \f[
      74             : D_{ij}^2 = \sum_{\alpha} (X^i_\alpha - X^j_\alpha)^2 = \sum_\alpha (X^i_\alpha)^2 + (X^j_\alpha)^2 - 2X^i_\alpha X^j_\alpha
      75             : \f]
      76             : 
      77             : We can use this expression and matrix algebra to calculate multiple distances at once.  For instance if we have three points,
      78             : \f$\mathbf{X}\f$, we can write distances between them as:
      79             : 
      80             : \f{eqnarray*}{
      81             : D^2(\mathbf{X}) &=& \left[ \begin{array}{ccc}
      82             : 0 & d_{12}^2 & d_{13}^2 \\
      83             : d_{12}^2 & 0 & d_{23}^2 \\
      84             : d_{13}^2 & d_{23}^2 & 0
      85             : \end{array}\right] \\
      86             : &=&
      87             : \sum_\alpha \left[ \begin{array}{ccc}
      88             : (X^1_\alpha)^2 & (X^1_\alpha)^2 & (X^1_\alpha)^2 \\
      89             : (X^2_\alpha)^2 & (X^2_\alpha)^2 & (X^2_\alpha)^2 \\
      90             : (X^3_\alpha)^2 & (X^3_\alpha)^2 & (X^3_\alpha)^2 \\
      91             : \end{array}\right]
      92             :  + \sum_\alpha \left[ \begin{array}{ccc}
      93             : (X^1_\alpha)^2 & (X^2_\alpha)^2 & (X^3_\alpha)^2 \\
      94             : (X^1_\alpha)^2 & (X^2_\alpha)^2 & (X^3_\alpha)^2 \\
      95             : (X^1_\alpha)^2 & (X^2_\alpha)^2 & (X^3_\alpha)^2 \\
      96             : \end{array}\right]
      97             : - 2 \sum_\alpha \left[ \begin{array}{ccc}
      98             : X^1_\alpha X^1_\alpha & X^1_\alpha X^2_\alpha & X^1_\alpha X^3_\alpha \\
      99             : X^2_\alpha X^1_\alpha & X^2_\alpha X^2_\alpha & X^2_\alpha X^3_\alpha \\
     100             : X^1_\alpha X^3_\alpha & X^3_\alpha X^2_\alpha & X^3_\alpha X^3_\alpha
     101             : \end{array}\right] \nonumber \\
     102             : &=& \mathbf{c 1^T} + \mathbf{1 c^T} - 2 \sum_\alpha \mathbf{x}_a \mathbf{x}^T_a =  \mathbf{c 1^T} + \mathbf{1 c^T} - 2\mathbf{X X^T}
     103             : \f}
     104             : 
     105             : This last equation can be extended to situations when we have more than three points.  In it \f$\mathbf{X}\f$ is a matrix that has
     106             : one high-dimensional point on each of its rows and \f$\mathbf{X^T}\f$ is its transpose.  \f$\mathbf{1}\f$ is an \f$M \times 1\f$ vector
     107             : of ones and \f$\mathbf{c}\f$ is a vector with components given by:
     108             : 
     109             : \f[
     110             : c_i = \sum_\alpha (x_\alpha^i)^2
     111             : \f]
     112             : 
     113             : These quantities are the diagonal elements of \f$\mathbf{X X^T}\f$, which is a dot product or Gram Matrix that contains the
     114             : dot product of the vector \f$X_i\f$ with the vector \f$X_j\f$ in element \f$i,j\f$.
     115             : 
     116             : In classical scaling we introduce a centering matrix \f$\mathbf{J}\f$ that is given by:
     117             : 
     118             : \f[
     119             : \mathbf{J} = \mathbf{I} - \frac{1}{M} \mathbf{11^T}
     120             : \f]
     121             : 
     122             : where \f$\mathbf{I}\f$ is the identity.  Multiplying the equations above from the front and back by this matrix and a factor of a \f$-\frac{1}{2}\f$ gives:
     123             : 
     124             : \f{eqnarray*}{
     125             :  -\frac{1}{2} \mathbf{J} \mathbf{D}^2(\mathbf{X}) \mathbf{J} &=& -\frac{1}{2}\mathbf{J}( \mathbf{c 1^T} + \mathbf{1 c^T} - 2\mathbf{X X^T})\mathbf{J} \\
     126             :  &=& -\frac{1}{2}\mathbf{J c 1^T J} - \frac{1}{2} \mathbf{J 1 c^T J} + \frac{1}{2} \mathbf{J}(2\mathbf{X X^T})\mathbf{J} \\
     127             :  &=& \mathbf{ J X X^T J } = \mathbf{X X^T } \label{eqn:scaling}
     128             : \f}
     129             : 
     130             : The fist two terms in this expression disappear because \f$\mathbf{1^T J}=\mathbf{J 1} =\mathbf{0}\f$, where \f$\mathbf{0}\f$
     131             : is a matrix containing all zeros.  In the final step meanwhile we use the fact that the matrix of squared distances will not
     132             : change when we translate all the points.  We can thus assume that the mean value, \f$\mu\f$, for each of the components, \f$\alpha\f$:
     133             : \f[
     134             : \mu_\alpha = \frac{1}{M} \sum_{i=1}^N \mathbf{X}^i_\alpha
     135             : \f]
     136             : is equal to 0 so the columns of \f$\mathbf{X}\f$ add up to 0.  This in turn means that each of the columns of
     137             : \f$\mathbf{X X^T}\f$ adds up to zero, which is what allows us to write \f$\mathbf{ J X X^T J } = \mathbf{X X^T }\f$.
     138             : 
     139             : The matrix of squared distances is symmetric and positive-definite we can thus use the spectral decomposition to decompose it as:
     140             : 
     141             : \f[
     142             : \Phi= \mathbf{V} \Lambda \mathbf{V}^T
     143             : \f]
     144             : 
     145             : Furthermore, because the matrix we are diagonalizing, \f$\mathbf{X X^T}\f$, is the product of a matrix and its transpose
     146             : we can use this decomposition to write:
     147             : 
     148             : \f[
     149             : \mathbf{X} =\mathbf{V} \Lambda^\frac{1}{2}
     150             : \f]
     151             : 
     152             : Much as in PCA there are generally a small number of large eigenvalues in \f$\Lambda\f$ and many small eigenvalues.
     153             : We can safely use only the large eigenvalues and their corresponding eigenvectors to express the relationship between
     154             : the coordinates \f$\mathbf{X}\f$.  This gives us our set of low-dimensional projections.
     155             : 
     156             : This derivation makes a number of assumptions about the how the low dimensional points should best be arranged to minimize
     157             : the stress. If you use an interactive optimization algorithm such as SMACOF you may thus be able to find a better
     158             : (lower-stress) projection of the points.  For more details on the assumptions made
     159             : see <a href="http://quest4rigor.com/tag/multidimensional-scaling/"> this website.</a>
     160             : */
     161             : //+ENDPLUMEDOC
     162             : 
     163             : namespace PLMD {
     164             : namespace dimred {
     165             : 
     166             : class ClassicalMultiDimensionalScaling : public ActionShortcut {
     167             : public:
     168             :   static void registerKeywords( Keywords& keys );
     169             :   explicit ClassicalMultiDimensionalScaling( const ActionOptions& ao );
     170             : };
     171             : 
     172             : PLUMED_REGISTER_ACTION(ClassicalMultiDimensionalScaling,"CLASSICAL_MDS")
     173             : 
     174          13 : void ClassicalMultiDimensionalScaling::registerKeywords( Keywords& keys ) {
     175          13 :   ActionShortcut::registerKeywords( keys );
     176          26 :   keys.add("compulsory","ARG","the arguments that you would like to make the histogram for");
     177          26 :   keys.add("compulsory","NLOW_DIM","number of low-dimensional coordinates required");
     178          26 :   keys.setValueDescription("matrix","the low dimensional projections for the input data points");
     179          52 :   keys.needsAction("TRANSPOSE"); keys.needsAction("DISSIMILARITIES"); keys.needsAction("MATRIX_VECTOR_PRODUCT"); keys.needsAction("VSTACK");
     180          52 :   keys.needsAction("SUM"); keys.needsAction("CUSTOM"); keys.needsAction("OUTER_PRODUCT"); keys.needsAction("DIAGONALIZE");
     181          13 : }
     182             : 
     183           7 : ClassicalMultiDimensionalScaling::ClassicalMultiDimensionalScaling( const ActionOptions& ao):
     184             :   Action(ao),
     185           7 :   ActionShortcut(ao)
     186             : {
     187             :   // Find the argument name
     188          14 :   std::string argn; parse("ARG",argn); std::string dissimilarities="";
     189           7 :   ActionShortcut* as = plumed.getActionSet().getShortcutActionWithLabel( argn );
     190           7 :   if( !as ) error("found no action with name " + argn );
     191           7 :   if( as->getName()!="COLLECT_FRAMES" ) {
     192           1 :     if( as->getName().find("LANDMARK_SELECT")==std::string::npos ) {
     193           0 :       error("found no COLLECT_FRAMES or LANDMARK_SELECT action with label " + argn );
     194             :     } else {
     195           1 :       ActionWithValue* dissims = plumed.getActionSet().selectWithLabel<ActionWithValue*>( argn + "_sqrdissims");
     196           2 :       if( dissims ) dissimilarities = argn + "_sqrdissims";
     197             :     }
     198             :   }
     199           7 :   if( dissimilarities.length()==0 ) {
     200           6 :     dissimilarities = getShortcutLabel() + "_mat";
     201             :     // Transpose matrix of stored data values
     202          12 :     readInputLine( argn + "_dataT: TRANSPOSE ARG=" + argn + "_data");
     203             :     // Calculate the dissimilarity matrix
     204          12 :     readInputLine( getShortcutLabel() + "_mat: DISSIMILARITIES SQUARED ARG=" + argn + "_data," + argn + "_dataT");
     205             :   }
     206             :   // Center the matrix
     207             :   // Step 1: calculate the sum of the rows and duplicate them into a matrix
     208          14 :   readInputLine( getShortcutLabel() + "_rsums: MATRIX_VECTOR_PRODUCT ARG=" + dissimilarities + "," + argn + "_ones" );
     209          14 :   readInputLine( getShortcutLabel() + "_nones: SUM ARG=" + argn + "_ones PERIODIC=NO");
     210          14 :   readInputLine( getShortcutLabel() + "_rsumsn: CUSTOM ARG=" + getShortcutLabel() + "_rsums," + getShortcutLabel() + "_nones FUNC=x/y PERIODIC=NO");
     211          14 :   readInputLine( getShortcutLabel() + "_rsummat: OUTER_PRODUCT ARG=" + getShortcutLabel() + "_rsumsn," + argn + "_ones");
     212             :   // Step 2: Multiply matrix by -0.5 and subtract row sums
     213          14 :   readInputLine( getShortcutLabel() + "_int: CUSTOM ARG=" + getShortcutLabel() + "_rsummat," + dissimilarities + " FUNC=-0.5*y+0.5*x PERIODIC=NO");
     214             :   // Step 3: Calculate column sums for new matrix and duplicate them into a matrix
     215          14 :   readInputLine( getShortcutLabel() + "_intT: TRANSPOSE ARG=" + getShortcutLabel() + "_int");
     216          14 :   readInputLine( getShortcutLabel() + "_csums: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_intT," + argn + "_ones" );
     217          14 :   readInputLine( getShortcutLabel() + "_csumsn: CUSTOM ARG=" + getShortcutLabel() + "_csums," + getShortcutLabel() + "_nones FUNC=x/y PERIODIC=NO");
     218          14 :   readInputLine( getShortcutLabel() + "_csummat: OUTER_PRODUCT ARG=" + getShortcutLabel() + "_csumsn," + argn + "_ones");
     219             :   // Step 4: subtract the column sums
     220          14 :   readInputLine( getShortcutLabel() + "_cmat: CUSTOM ARG=" + getShortcutLabel() + "_csummat," + getShortcutLabel() + "_intT FUNC=y-x PERIODIC=NO");
     221             :   // And generate the multidimensional scaling projection
     222           7 :   unsigned ndim; parse("NLOW_DIM",ndim);
     223          19 :   std::string vecstr="1"; for(unsigned i=1; i<ndim; ++i) { std::string num; Tools::convert( i+1, num ); vecstr += "," + num; }
     224          14 :   readInputLine( getShortcutLabel() + "_eig: DIAGONALIZE ARG=" + getShortcutLabel() + "_cmat VECTORS=" + vecstr );
     225          20 :   for(unsigned i=0; i<ndim; ++i) {
     226          13 :     std::string num; Tools::convert( i+1, num );
     227          26 :     readInputLine( getShortcutLabel() + "-" +  num + ": CUSTOM ARG=" + getShortcutLabel() + "_eig.vals-" + num + "," + getShortcutLabel() + "_eig.vecs-" + num + " FUNC=sqrt(x)*y PERIODIC=NO");
     228             :   }
     229           7 :   std::string eigvec_args = " ARG=" + getShortcutLabel() + "-1";
     230             :   // The final output is a stack of all the low dimensional coordinates
     231          19 :   for(unsigned i=1; i<ndim; ++i) { std::string num; Tools::convert( i+1, num ); eigvec_args += "," + getShortcutLabel() + "-" + num; }
     232          14 :   readInputLine( getShortcutLabel() + ": VSTACK" + eigvec_args );
     233           7 : }
     234             : 
     235             : }
     236             : }