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 [isometry](http://en.wikipedia.org/wiki/Isometry)
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 $M$ $D$-dimensional points, $\mathbf{X}$,
40 : and we can calculate dissimilarities between pairs them, $D_{ij}$, we can, with an MDS calculation, try to create $M$ projections,
41 : $\mathbf{x}$, of the high dimensionality points in a $d$-dimensional linear space by trying to arrange the projections so that the
42 : Euclidean distances between pairs of them, $d_{ij}$, resemble the dissimilarities between the high dimensional points. In short we minimize:
43 :
44 : $$
45 : \chi^2 = \sum_{i \ne j} \left( D_{ij} - d_{ij} \right)^2
46 : $$
47 :
48 : where $D_{ij}$ is the distance between point $X^{i}$ and point $X^{j}$ and $d_{ij}$ is the distance between the projection
49 : of $X^{i}$, $x^i$, and the projection of $X^{j}$, $x^j$. A tutorial on this approach can be used to analyze simulations
50 : can be found in [this tutorial](https://www.plumed-tutorials.org/lessons/21/006/data/DIMENSIONALITY.html) and in the following [short video](https://www.youtube.com/watch?v=ofC2qz0_9_A&feature=youtu.be).
51 :
52 : ## 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 : ```plumed
58 : ff: COLLECT_FRAMES ATOMS=1-256
59 : mds: CLASSICAL_MDS ARG=ff NLOW_DIM=2
60 : weights: CUSTOM ARG=ff.logweights FUNC=exp(x) PERIODIC=NO
61 : DUMPPDB ATOM_INDICES=1-256 ATOMS=ff_data ARG=mds,weights FILE=embed.pdb
62 : ```
63 :
64 : By contrast the following input instructs PLUMED to calculate the distances between atoms by taking the differences in five torsional angles.
65 : The MDS algorithm is then used to arrange a set of points in a low dimensional space in a way that reproduces these differences.
66 :
67 : ```plumed
68 : phi1: TORSION ATOMS=1,2,3,4
69 : phi2: TORSION ATOMS=5,6,7,8
70 : phi3: TORSION ATOMS=9,10,11,12
71 : phi4: TORSION ATOMS=13,14,15,16
72 : phi5: TORSION ATOMS=17,18,19,20
73 :
74 : angles: COLLECT_FRAMES ARG=phi1,phi2,phi3,phi4,phi5
75 : mds: CLASSICAL_MDS ARG=angles NLOW_DIM=2
76 : weights: CUSTOM ARG=angles.logweights FUNC=exp(x) PERIODIC=NO
77 : DUMPVECTOR ARG=mds,weights FILE=list_embed
78 : ```
79 :
80 : The following section is for people who are interested in how this method works in detail. A solid understanding of this material is
81 : not necessary to use MDS.
82 :
83 : ## Method of optimization
84 :
85 : The stress function can be minimized using the conjugate gradients or steepest descent optimization algorithms that are implemented in [ARRANGE_POINTS](ARRANGE_POINTS.md).
86 : However, it is more common to do this minimization using a technique known as classical scaling. Classical scaling works by
87 : recognizing that each of the distances $D_{ij}$ in the above sum can be written as:
88 :
89 : $$
90 : 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
91 : $$
92 :
93 : We can use this expression and matrix algebra to calculate multiple distances at once. For instance if we have three points,
94 : $\mathbf{X}$, we can write distances between them as:
95 :
96 : $$
97 : \begin{aligned}
98 : D^2(\mathbf{X}) &= \left[ \begin{array}{ccc}
99 : 0 & d_{12}^2 & d_{13}^2 \\
100 : d_{12}^2 & 0 & d_{23}^2 \\
101 : d_{13}^2 & d_{23}^2 & 0
102 : \end{array}\right] \\
103 : &=
104 : \sum_\alpha \left[ \begin{array}{ccc}
105 : (X^1_\alpha)^2 & (X^1_\alpha)^2 & (X^1_\alpha)^2 \\
106 : (X^2_\alpha)^2 & (X^2_\alpha)^2 & (X^2_\alpha)^2 \\
107 : (X^3_\alpha)^2 & (X^3_\alpha)^2 & (X^3_\alpha)^2 \\
108 : \end{array}\right]
109 : + \sum_\alpha \left[ \begin{array}{ccc}
110 : (X^1_\alpha)^2 & (X^2_\alpha)^2 & (X^3_\alpha)^2 \\
111 : (X^1_\alpha)^2 & (X^2_\alpha)^2 & (X^3_\alpha)^2 \\
112 : (X^1_\alpha)^2 & (X^2_\alpha)^2 & (X^3_\alpha)^2 \\
113 : \end{array}\right]
114 : - 2 \sum_\alpha \left[ \begin{array}{ccc}
115 : X^1_\alpha X^1_\alpha & X^1_\alpha X^2_\alpha & X^1_\alpha X^3_\alpha \\
116 : X^2_\alpha X^1_\alpha & X^2_\alpha X^2_\alpha & X^2_\alpha X^3_\alpha \\
117 : X^1_\alpha X^3_\alpha & X^3_\alpha X^2_\alpha & X^3_\alpha X^3_\alpha
118 : \end{array}\right] \nonumber \\
119 : &= \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}
120 : \end{aligned}
121 : $$
122 :
123 : This last equation can be extended to situations when we have more than three points. In it $\mathbf{X}$ is a matrix that has
124 : one high-dimensional point on each of its rows and $\mathbf{X^T}$ is its transpose. $\mathbf{1}$ is an $M \times 1$ vector
125 : of ones and $\mathbf{c}$ is a vector with components given by:
126 :
127 : $$
128 : c_i = \sum_\alpha (x_\alpha^i)^2
129 : $$
130 :
131 : These quantities are the diagonal elements of $\mathbf{X X^T}$, which is a dot product or Gram Matrix that contains the
132 : dot product of the vector $X_i$ with the vector $X_j$ in element $i,j$.
133 :
134 : In classical scaling we introduce a centering matrix $\mathbf{J}$ that is given by:
135 :
136 : $$
137 : \mathbf{J} = \mathbf{I} - \frac{1}{M} \mathbf{11^T}
138 : $$
139 :
140 : where $\mathbf{I}$ is the identity. Multiplying the equations above from the front and back by this matrix and a factor of a $-\frac{1}{2}$ gives:
141 :
142 : $$
143 : \begin{aligned}
144 : -\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} \\
145 : &= -\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} \\
146 : &= \mathbf{ J X X^T J } = \mathbf{X X^T } \label{eqn:scaling}
147 : \end{aligned}
148 : $$
149 :
150 : The fist two terms in this expression disappear because $\mathbf{1^T J}=\mathbf{J 1} =\mathbf{0}$, where $\mathbf{0}$
151 : is a matrix containing all zeros. In the final step meanwhile we use the fact that the matrix of squared distances will not
152 : change when we translate all the points. We can thus assume that the mean value, $\mu$, for each of the components, $\alpha$:
153 :
154 : $$
155 : \mu_\alpha = \frac{1}{M} \sum_{i=1}^N \mathbf{X}^i_\alpha
156 : $$
157 :
158 : is equal to 0 so the columns of $\mathbf{X}$ add up to 0. This in turn means that each of the columns of
159 : $\mathbf{X X^T}$ adds up to zero, which is what allows us to write $\mathbf{ J X X^T J } = \mathbf{X X^T }$.
160 :
161 : The matrix of squared distances is symmetric and positive-definite we can thus use the spectral decomposition to decompose it as:
162 :
163 : $$
164 : \Phi= \mathbf{V} \Lambda \mathbf{V}^T
165 : $$
166 :
167 : Furthermore, because the matrix we are diagonalizing, $\mathbf{X X^T}$, is the product of a matrix and its transpose
168 : we can use this decomposition to write:
169 :
170 : $$
171 : \mathbf{X} =\mathbf{V} \Lambda^\frac{1}{2}
172 : $$
173 :
174 : Much as in PCA there are generally a small number of large eigenvalues in $\Lambda$ and many small eigenvalues.
175 : We can safely use only the large eigenvalues and their corresponding eigenvectors to express the relationship between
176 : the coordinates $\mathbf{X}$. This gives us our set of low-dimensional projections.
177 :
178 : This derivation makes a number of assumptions about the how the low dimensional points should best be arranged to minimize
179 : the stress. If you use an interactive optimization algorithm such as SMACOF you may thus be able to find a better
180 : (lower-stress) projection of the points. For more details on the assumptions made
181 : see [this website](http://quest4rigor.com/tag/multidimensional-scaling/).
182 : */
183 : //+ENDPLUMEDOC
184 :
185 : namespace PLMD {
186 : namespace dimred {
187 :
188 : class ClassicalMultiDimensionalScaling : public ActionShortcut {
189 : public:
190 : static void registerKeywords( Keywords& keys );
191 : explicit ClassicalMultiDimensionalScaling( const ActionOptions& ao );
192 : };
193 :
194 : PLUMED_REGISTER_ACTION(ClassicalMultiDimensionalScaling,"CLASSICAL_MDS")
195 :
196 13 : void ClassicalMultiDimensionalScaling::registerKeywords( Keywords& keys ) {
197 13 : ActionShortcut::registerKeywords( keys );
198 13 : keys.add("compulsory","ARG","the arguments that you would like to make the histogram for");
199 13 : keys.add("compulsory","NLOW_DIM","number of low-dimensional coordinates required");
200 26 : keys.setValueDescription("matrix","the low dimensional projections for the input data points");
201 13 : keys.needsAction("TRANSPOSE");
202 13 : keys.needsAction("DISSIMILARITIES");
203 13 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
204 13 : keys.needsAction("VSTACK");
205 13 : keys.needsAction("SUM");
206 13 : keys.needsAction("CUSTOM");
207 13 : keys.needsAction("OUTER_PRODUCT");
208 13 : keys.needsAction("DIAGONALIZE");
209 13 : }
210 :
211 7 : ClassicalMultiDimensionalScaling::ClassicalMultiDimensionalScaling( const ActionOptions& ao):
212 : Action(ao),
213 7 : ActionShortcut(ao) {
214 : // Find the argument name
215 : std::string argn;
216 7 : parse("ARG",argn);
217 7 : std::string dissimilarities="";
218 7 : ActionShortcut* as = plumed.getActionSet().getShortcutActionWithLabel( argn );
219 7 : if( !as ) {
220 0 : error("found no action with name " + argn );
221 : }
222 7 : if( as->getName()!="COLLECT_FRAMES" ) {
223 1 : if( as->getName().find("LANDMARK_SELECT")==std::string::npos ) {
224 0 : error("found no COLLECT_FRAMES or LANDMARK_SELECT action with label " + argn );
225 : } else {
226 1 : ActionWithValue* dissims = plumed.getActionSet().selectWithLabel<ActionWithValue*>( argn + "_sqrdissims");
227 1 : if( dissims ) {
228 2 : dissimilarities = argn + "_sqrdissims";
229 : }
230 : }
231 : }
232 7 : if( dissimilarities.length()==0 ) {
233 6 : dissimilarities = getShortcutLabel() + "_mat";
234 : // Transpose matrix of stored data values
235 12 : readInputLine( argn + "_dataT: TRANSPOSE ARG=" + argn + "_data");
236 : // Calculate the dissimilarity matrix
237 12 : readInputLine( getShortcutLabel() + "_mat: DISSIMILARITIES SQUARED ARG=" + argn + "_data," + argn + "_dataT");
238 : }
239 : // Center the matrix
240 : // Step 1: calculate the sum of the rows and duplicate them into a matrix
241 14 : readInputLine( getShortcutLabel() + "_rsums: MATRIX_VECTOR_PRODUCT ARG=" + dissimilarities + "," + argn + "_ones" );
242 14 : readInputLine( getShortcutLabel() + "_nones: SUM ARG=" + argn + "_ones PERIODIC=NO");
243 14 : readInputLine( getShortcutLabel() + "_rsumsn: CUSTOM ARG=" + getShortcutLabel() + "_rsums," + getShortcutLabel() + "_nones FUNC=x/y PERIODIC=NO");
244 14 : readInputLine( getShortcutLabel() + "_rsummat: OUTER_PRODUCT ARG=" + getShortcutLabel() + "_rsumsn," + argn + "_ones");
245 : // Step 2: Multiply matrix by -0.5 and subtract row sums
246 14 : readInputLine( getShortcutLabel() + "_int: CUSTOM ARG=" + getShortcutLabel() + "_rsummat," + dissimilarities + " FUNC=-0.5*y+0.5*x PERIODIC=NO");
247 : // Step 3: Calculate column sums for new matrix and duplicate them into a matrix
248 14 : readInputLine( getShortcutLabel() + "_intT: TRANSPOSE ARG=" + getShortcutLabel() + "_int");
249 14 : readInputLine( getShortcutLabel() + "_csums: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_intT," + argn + "_ones" );
250 14 : readInputLine( getShortcutLabel() + "_csumsn: CUSTOM ARG=" + getShortcutLabel() + "_csums," + getShortcutLabel() + "_nones FUNC=x/y PERIODIC=NO");
251 14 : readInputLine( getShortcutLabel() + "_csummat: OUTER_PRODUCT ARG=" + getShortcutLabel() + "_csumsn," + argn + "_ones");
252 : // Step 4: subtract the column sums
253 14 : readInputLine( getShortcutLabel() + "_cmat: CUSTOM ARG=" + getShortcutLabel() + "_csummat," + getShortcutLabel() + "_intT FUNC=y-x PERIODIC=NO");
254 : // And generate the multidimensional scaling projection
255 : unsigned ndim;
256 7 : parse("NLOW_DIM",ndim);
257 7 : std::string vecstr="1";
258 13 : for(unsigned i=1; i<ndim; ++i) {
259 : std::string num;
260 6 : Tools::convert( i+1, num );
261 12 : vecstr += "," + num;
262 : }
263 14 : readInputLine( getShortcutLabel() + "_eig: DIAGONALIZE ARG=" + getShortcutLabel() + "_cmat VECTORS=" + vecstr );
264 20 : for(unsigned i=0; i<ndim; ++i) {
265 : std::string num;
266 13 : Tools::convert( i+1, num );
267 26 : readInputLine( getShortcutLabel() + "-" + num + ": CUSTOM ARG=" + getShortcutLabel() + "_eig.vecs-" + num + "," + getShortcutLabel() + "_eig.vals-" + num + " FUNC=sqrt(y)*x PERIODIC=NO");
268 : }
269 7 : std::string eigvec_args = " ARG=" + getShortcutLabel() + "-1";
270 : // The final output is a stack of all the low dimensional coordinates
271 13 : for(unsigned i=1; i<ndim; ++i) {
272 : std::string num;
273 6 : Tools::convert( i+1, num );
274 12 : eigvec_args += "," + getShortcutLabel() + "-" + num;
275 : }
276 14 : readInputLine( getShortcutLabel() + ": VSTACK" + eigvec_args );
277 7 : }
278 :
279 : }
280 : }
|
