diff --git a/.travis.yml b/.travis.yml
index e7429211ce476211fb3b10f24d62b0e50fdb7339..729c7430c17ac025898bf1632bf70694dc7f2451 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -66,7 +66,7 @@ install:
- if test -n "$LD_LIBRARY_PATH" ; then PLUMED_LDFLAGS="-Wl,-rpath,$LD_LIBRARY_PATH" ; fi
script:
# we set all the optional modules on
- - touch src/crystallization.on src/manyrestraints.on
+ - touch src/crystallization.on src/manyrestraints.on src/dimred.on
# BUILD:
# this is done only if PLUMED_CXX is defined
diff --git a/regtest/dimred/rt-mds/Makefile b/regtest/dimred/rt-mds/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..3703b27cea227aa053fb6d1d73f861e4384dbcee
--- /dev/null
+++ b/regtest/dimred/rt-mds/Makefile
@@ -0,0 +1 @@
+include ../../scripts/test.make
diff --git a/regtest/dimred/rt-mds/analysis.0.list_embed.reference b/regtest/dimred/rt-mds/analysis.0.list_embed.reference
new file mode 100644
index 0000000000000000000000000000000000000000..ce29c25f6a5ecf50ae070443c00face1089a9d06
--- /dev/null
+++ b/regtest/dimred/rt-mds/analysis.0.list_embed.reference
@@ -0,0 +1,101 @@
+#! FIELDS mds.1 mds.2
+ 0.0873 0.0013
+ 0.0535 0.0044
+ 0.0150 0.0044
+ -0.0136 0.0019
+ -0.0308 -0.0012
+ -0.0355 -0.0018
+ -0.0090 -0.0015
+ 0.0145 0.0004
+ 0.0191 0.0022
+ 0.0278 0.0015
+ 0.0091 -0.0004
+ -0.0359 -0.0012
+ -0.0540 0.0004
+ -0.0322 -0.0034
+ -0.0205 -0.0150
+ -0.0282 -0.0289
+ -0.0530 -0.0358
+ -0.0769 -0.0312
+ -0.0893 -0.0259
+ -0.0697 -0.0246
+ -0.0301 -0.0196
+ 0.0093 -0.0072
+ 0.0274 0.0009
+ 0.0301 0.0016
+ 0.0265 -0.0033
+ 0.0433 -0.0053
+ 0.0514 -0.0019
+ 0.0222 0.0031
+ -0.0129 0.0061
+ -0.0322 0.0084
+ -0.0355 0.0080
+ -0.0337 0.0067
+ -0.0308 0.0047
+ -0.0253 0.0016
+ -0.0132 -0.0021
+ 0.0040 -0.0031
+ 0.0107 0.0011
+ 0.0299 0.0060
+ 0.0540 0.0070
+ 0.0650 0.0039
+ 0.0415 0.0025
+ 0.0074 0.0033
+ -0.0185 0.0059
+ -0.0413 0.0073
+ -0.0548 0.0074
+ -0.0385 0.0041
+ -0.0383 -0.0006
+ -0.0584 -0.0008
+ -0.0584 0.0023
+ -0.0553 0.0090
+ -0.0772 0.0146
+ -0.0757 0.0157
+ -0.0644 0.0144
+ -0.0816 0.0141
+ -0.1045 0.0148
+ -0.0867 0.0111
+ -0.0365 0.0014
+ 0.0113 -0.0067
+ 0.0272 -0.0074
+ 0.0381 -0.0012
+ 0.0548 0.0021
+ 0.0470 0.0012
+ 0.0227 -0.0019
+ 0.0132 -0.0022
+ 0.0040 0.0004
+ -0.0045 0.0008
+ -0.0145 0.0006
+ -0.0116 0.0025
+ -0.0031 0.0044
+ -0.0148 0.0086
+ -0.0318 0.0109
+ -0.0166 0.0078
+ 0.0154 0.0036
+ 0.0451 0.0018
+ 0.0364 0.0021
+ 0.0125 0.0036
+ 0.0330 0.0045
+ 0.0605 0.0048
+ 0.0417 0.0041
+ 0.0115 0.0011
+ 0.0257 -0.0018
+ 0.0451 -0.0029
+ 0.0387 -0.0032
+ 0.0254 -0.0017
+ 0.0311 -0.0001
+ 0.0448 0.0006
+ 0.0499 0.0020
+ 0.0543 0.0037
+ 0.0486 0.0054
+ 0.0459 0.0051
+ 0.0715 -0.0001
+ 0.0764 -0.0063
+ 0.0603 -0.0119
+ 0.0378 -0.0153
+ 0.0214 -0.0124
+ 0.0063 -0.0091
+ -0.0095 -0.0016
+ -0.0259 0.0036
+ -0.0234 0.0061
+ 0.0012 0.0060
diff --git a/regtest/dimred/rt-mds/config b/regtest/dimred/rt-mds/config
new file mode 100644
index 0000000000000000000000000000000000000000..57d886c44c9b5b615e6e31bf117a02928408814b
--- /dev/null
+++ b/regtest/dimred/rt-mds/config
@@ -0,0 +1,2 @@
+type=simplemd
+
diff --git a/regtest/dimred/rt-mds/in b/regtest/dimred/rt-mds/in
new file mode 100644
index 0000000000000000000000000000000000000000..030417936fdc6c7f1f4d0edff090893b1485f13d
--- /dev/null
+++ b/regtest/dimred/rt-mds/in
@@ -0,0 +1,11 @@
+inputfile input.xyz
+outputfile output.xyz
+temperature 0.2
+tstep 0.005
+friction 1
+forcecutoff 2.5
+listcutoff 3.0
+ndim 2
+nstep 2000
+nconfig 1000 trajectory.xyz
+nstat 1000 energies.dat
diff --git a/regtest/dimred/rt-mds/input.xyz b/regtest/dimred/rt-mds/input.xyz
new file mode 100644
index 0000000000000000000000000000000000000000..9e1e59174c1ce7a8d31fa8220d64ea05a9edbc14
--- /dev/null
+++ b/regtest/dimred/rt-mds/input.xyz
@@ -0,0 +1,9 @@
+7
+100. 100. 100.
+Ar 7.3933470660 -2.6986483924 0.0000000000
+Ar 7.8226765198 -0.7390907295 0.0000000000
+Ar 7.1014969839 -1.6164766614 0.0000000000
+Ar 8.2357184242 -1.7097824975 0.0000000000
+Ar 6.7372520842 -0.5111536183 0.0000000000
+Ar 6.3777119489 -2.4640437401 0.0000000000
+Ar 5.9900631495 -1.3385375043 0.0000000000
diff --git a/regtest/dimred/rt-mds/list_embed.reference b/regtest/dimred/rt-mds/list_embed.reference
new file mode 100644
index 0000000000000000000000000000000000000000..2f4fb6096a9baef0d55b0c24a64cc504a1b9021a
--- /dev/null
+++ b/regtest/dimred/rt-mds/list_embed.reference
@@ -0,0 +1,101 @@
+#! FIELDS mds.1 mds.2
+ 0.0594 0.0053
+ 0.0683 0.0034
+ 0.0705 0.0006
+ 0.0729 -0.0019
+ 0.0874 -0.0069
+ 0.0961 -0.0095
+ 0.0867 -0.0062
+ 0.0719 -0.0008
+ 0.0575 0.0031
+ 0.0552 0.0034
+ 0.0682 -0.0002
+ 0.0781 -0.0024
+ 0.0785 -0.0025
+ 0.0763 -0.0027
+ 0.0755 -0.0041
+ 0.0792 -0.0069
+ 0.0778 -0.0062
+ 0.0727 -0.0025
+ 0.0726 0.0014
+ 0.0756 0.0030
+ 0.0821 0.0017
+ 0.0816 0.0001
+ 0.0774 0.0008
+ 0.0766 0.0008
+ 0.0838 -0.0005
+ 0.0931 -0.0016
+ 0.0936 -0.0010
+ 0.0923 -0.0009
+ 0.0906 -0.0008
+ 0.0706 0.0025
+ 0.0664 0.0050
+ 0.0901 0.0019
+ 0.1059 -0.0015
+ 0.0917 0.0006
+ 0.0876 0.0021
+ 0.1008 -0.0001
+ 0.0955 -0.0014
+ 0.0858 -0.0018
+ 0.0894 -0.0005
+ 0.0922 0.0014
+ 0.0917 0.0010
+ 0.0925 -0.0043
+ 0.0967 -0.0102
+ 0.1073 -0.0096
+ 0.1065 -0.0047
+ 0.0935 -0.0003
+ 0.0668 0.0049
+ 0.0411 0.0087
+ 0.0264 0.0110
+ 0.0254 0.0125
+ 0.0503 0.0085
+ 0.0855 0.0007
+ 0.0975 -0.0028
+ 0.0983 -0.0035
+ 0.0997 -0.0044
+ 0.0877 -0.0021
+ 0.0501 0.0053
+ 0.0104 0.0140
+ -0.0199 0.0221
+ -0.0416 0.0277
+ -0.0688 0.0338
+ -0.0858 0.0354
+ -0.0897 0.0313
+ -0.0920 0.0254
+ -0.1168 0.0227
+ -0.1447 0.0227
+ -0.1636 0.0221
+ -0.1461 0.0182
+ -0.1212 0.0139
+ -0.0775 0.0056
+ -0.0008 -0.0041
+ 0.0469 -0.0038
+ 0.0771 -0.0004
+ 0.0917 0.0005
+ 0.0656 0.0042
+ -0.0034 0.0115
+ -0.0901 0.0141
+ -0.1784 0.0069
+ -0.2253 -0.0077
+ -0.2372 -0.0248
+ -0.2415 -0.0448
+ -0.2562 -0.0539
+ -0.2901 -0.0433
+ -0.3173 -0.0215
+ -0.3117 0.0096
+ -0.2787 0.0317
+ -0.2319 0.0372
+ -0.1897 0.0305
+ -0.1592 0.0142
+ -0.1211 -0.0049
+ -0.0960 -0.0150
+ -0.0842 -0.0162
+ -0.0813 -0.0187
+ -0.0823 -0.0193
+ -0.0773 -0.0179
+ -0.0621 -0.0202
+ -0.0449 -0.0215
+ -0.0269 -0.0307
+ -0.0075 -0.0368
+ 0.0266 -0.0342
diff --git a/regtest/dimred/rt-mds/plumed.dat b/regtest/dimred/rt-mds/plumed.dat
new file mode 100755
index 0000000000000000000000000000000000000000..61954ad8fbb22eb59a5e6b575de00d96c8824b54
--- /dev/null
+++ b/regtest/dimred/rt-mds/plumed.dat
@@ -0,0 +1,29 @@
+UNITS NATURAL
+COM ATOMS=1-7 LABEL=com
+DISTANCE ATOMS=1,com LABEL=d1
+UPPER_WALLS ARG=d1 AT=2.0 KAPPA=100.
+DISTANCE ATOMS=2,com LABEL=d2
+UPPER_WALLS ARG=d2 AT=2.0 KAPPA=100.
+DISTANCE ATOMS=3,com LABEL=d3
+UPPER_WALLS ARG=d3 AT=2.0 KAPPA=100.
+DISTANCE ATOMS=4,com LABEL=d4
+UPPER_WALLS ARG=d4 AT=2.0 KAPPA=100.
+DISTANCE ATOMS=5,com LABEL=d5
+UPPER_WALLS ARG=d5 AT=2.0 KAPPA=100.
+DISTANCE ATOMS=6,com LABEL=d6
+UPPER_WALLS ARG=d6 AT=2.0 KAPPA=100.
+DISTANCE ATOMS=7,com LABEL=d7
+UPPER_WALLS ARG=d7 AT=2.0 KAPPA=100.
+
+COORDINATIONNUMBER SPECIES=1-7 MOMENTS=2-3 SWITCH={RATIONAL R_0=1.5 NN=8 MM=16} LABEL=c1
+
+oo: EUCLIDEAN_DISSIMILARITIES ARG=c1.moment-2,c1.moment-3 STRIDE=10 RUN=1000 REWEIGHT_TEMP=0.1 TEMP=0.2
+
+CLASSICAL_MDS ...
+ USE_OUTPUT_DATA_FROM=oo
+ NLOW_DIM=2
+ LABEL=mds
+... CLASSICAL_MDS
+
+OUTPUT_ANALYSIS_DATA_TO_COLVAR USE_OUTPUT_DATA_FROM=mds FILE=list_embed FMT=%8.4f
+# OUTPUT_ANALYSIS_DATA_TO_PDB USE_OUTPUT_DATA_FROM=mds FILE=embed FMT=%8.4f
diff --git a/src/.gitignore b/src/.gitignore
index 838ba6a71f6df12856a06c1dcd9cb7ca4f03f657..27396db47dd7ec10580bfda5483abc8e2b2b7097 100644
--- a/src/.gitignore
+++ b/src/.gitignore
@@ -11,6 +11,7 @@
!/config
!/core
!/crystallization
+!/dimred
!/function
!/generic
!/header.sh
diff --git a/src/dimred/.gitignore b/src/dimred/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..4888552338441b5d99a97ec7a2cb2ce448e9e4e6
--- /dev/null
+++ b/src/dimred/.gitignore
@@ -0,0 +1,9 @@
+/*
+# in this directory, only accept source, Makefile and README
+!/.gitignore
+!/*.c
+!/*.cpp
+!/*.h
+!/Makefile
+!/README
+!/module.type
diff --git a/src/dimred/ClassicalMultiDimensionalScaling.cpp b/src/dimred/ClassicalMultiDimensionalScaling.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9416579d8b2fa4ff3cf35cfb8aa9c4bd5fe4ca33
--- /dev/null
+++ b/src/dimred/ClassicalMultiDimensionalScaling.cpp
@@ -0,0 +1,213 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2012-2014 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "DimensionalityReductionBase.h"
+#include "core/ActionRegister.h"
+
+//+PLUMEDOC ANALYSIS CLASSICAL_MDS
+/*
+Create a low-dimensional projection of a trajectory using the classical multidimensional
+scaling algorithm.
+
+Multidimensional scaling (MDS) is similar to what is done when you make a map. You start with distances
+between London, Belfast, Paris and Dublin and then you try to arrange points on a piece of paper so that the (suitably scaled)
+distances between the points in your map representing each of those cities are related to the true distances between the cities.
+Stating this more mathematically MDS endeavors to find an <a href="http://en.wikipedia.org/wiki/Isometry">isometry</a>
+between points distributed in a high-dimensional space and a set of points distributed in a low-dimensional plane.
+In other words, if we have \f$M\f$ \f$D\f$-dimensional points, \f$\mathbf{X}\f$,
+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,
+\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
+Euclidean distances between pairs of them, \f$d_{ij}\f$, resemble the dissimilarities between the high dimensional points. In short we minimize:
+
+\f[
+\chi^2 = \sum_{i \ne j} \left( D_{ij} - d_{ij} \right)^2
+\f]
+
+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
+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 analyse simulations
+can be found in the tutorial \ref belfast-3 and in the following <a href="https://www.youtube.com/watch?v=ofC2qz0_9_A&feature=youtu.be" > short video.</a>
+
+\par Examples
+
+The following command instructs plumed to construct a classical multidimensional scaling projection of a trajectory.
+The RMSD distance between atoms 1-256 have moved is used to measure the distances in the high-dimensional space.
+
+\verbatim
+CLASSICAL_MDS ...
+ ATOMS=1-256
+ METRIC=OPTIMAL-FAST
+ USE_ALL_DATA
+ NLOW_DIM=2
+ OUTPUT_FILE=rmsd-embed
+... CLASSICAL_MDS
+\endverbatim
+
+The following section is for people who are interested in how this method works in detail. A solid understanding of this material is
+not necessary to use MDS.
+
+\section dim-sec Method of optimisation
+
+The stress function can be minimized using a standard optimization algorithm such as conjugate gradients or steepest descent.
+However, it is more common to do this minimization using a technique known as classical scaling. Classical scaling works by
+recognizing that each of the distances $D_{ij}$ in the above sum can be written as:
+
+\f[
+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
+\f]
+
+We can use this expression and matrix algebra to calculate multiple distances at once. For instance if we have three points,
+\f$\mathbf{X}\f$, we can write distances between them as:
+
+\f{eqnarray*}{
+D^2(\mathbf{X}) &=& \left[ \begin{array}{ccc}
+0 & d_{12}^2 & d_{13}^2 \\
+d_{12}^2 & 0 & d_{23}^2 \\
+d_{13}^2 & d_{23}^2 & 0
+\end{array}\right] \\
+&=&
+\sum_\alpha \left[ \begin{array}{ccc}
+(X^1_\alpha)^2 & (X^1_\alpha)^2 & (X^1_\alpha)^2 \\
+(X^2_\alpha)^2 & (X^2_\alpha)^2 & (X^2_\alpha)^2 \\
+(X^3_\alpha)^2 & (X^3_\alpha)^2 & (X^3_\alpha)^2 \\
+\end{array}\right]
+ + \sum_\alpha \left[ \begin{array}{ccc}
+(X^1_\alpha)^2 & (X^2_\alpha)^2 & (X^3_\alpha)^2 \\
+(X^1_\alpha)^2 & (X^2_\alpha)^2 & (X^3_\alpha)^2 \\
+(X^1_\alpha)^2 & (X^2_\alpha)^2 & (X^3_\alpha)^2 \\
+\end{array}\right]
+- 2 \sum_\alpha \left[ \begin{array}{ccc}
+X^1_\alpha X^1_\alpha & X^1_\alpha X^2_\alpha & X^1_\alpha X^3_\alpha \\
+X^2_\alpha X^1_\alpha & X^2_\alpha X^2_\alpha & X^2_\alpha X^3_\alpha \\
+X^1_\alpha X^3_\alpha & X^3_\alpha X^2_\alpha & X^3_\alpha X^3_\alpha
+\end{array}\right] \nonumber \\
+&=& \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}
+\f}
+
+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
+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
+of ones and \f$\mathbf{c}\f$ is a vector with components given by:
+
+\f[
+c_i = \sum_\alpha (x_\alpha^i)^2
+\f]
+
+These quantities are the diagonal elements of \f$\mathbf{X X^T}\f$, which is a dot product or Gram Matrix that contains the
+dot product of the vector \f$X_i\f$ with the vector \f$X_j\f$ in element \f$i,j\f$.
+
+In classical scaling we introduce a centering matrix \f$\mathbf{J}\f$ that is given by:
+
+\f[
+\mathbf{J} = \mathbf{I} - \frac{1}{M} \mathbf{11^T}
+\f]
+
+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:
+
+\f{eqnarray*}{
+ -\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} \\
+ &=& -\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} \\
+ &=& \mathbf{ J X X^T J } = \mathbf{X X^T } \label{eqn:scaling}
+\f}
+
+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$
+is a matrix containing all zeros. In the final step meanwhile we use the fact that the matrix of squared distances will not
+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$:
+\f[
+\mu_\alpha = \frac{1}{M} \sum_{i=1}^N \mathbf{X}^i_\alpha
+\f]
+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
+\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$.
+
+The matrix of squared distances is symmetric and positive-definite we can thus use the spectral decomposition to decompose it as:
+
+\f[
+\Phi= \mathbf{V} \Lambda \mathbf{V}^T
+\f]
+
+Furthermore, because the matrix we are diagonalizing, \f$\mathbf{X X^T}\f$, is the product of a matrix and its transpose
+we can use this decomposition to write:
+
+\f[
+\mathbf{X} =\mathbf{V} \Lambda^\frac{1}{2}
+\f]
+
+Much as in PCA there are generally a small number of large eigenvalues in \f$\Lambda\f$ and many small eigenvalues.
+We can safely use only the large eigenvalues and their corresponding eigenvectors to express the relationship between
+the coordinates \f$\mathbf{X}\f$. This gives us our set of low-dimensional projections.
+
+This derivation makes a number of assumptions about the how the low dimensional points should best be arranged to minimise
+the stress. If you use an interative optimization algorithm such as SMACOF you may thus be able to find a better
+(lower-stress) projection of the points. For more details on the assumptions made
+see <a href="http://quest4rigor.com/tag/multidimensional-scaling/"> this website.</a>
+*/
+//+ENDPLUMEDOC
+
+namespace PLMD {
+namespace dimred {
+
+class ClassicalMultiDimensionalScaling : public DimensionalityReductionBase {
+public:
+ static void registerKeywords( Keywords& keys );
+ ClassicalMultiDimensionalScaling( const ActionOptions& );
+ void calculateProjections( const Matrix<double>& , Matrix<double>& );
+};
+
+PLUMED_REGISTER_ACTION(ClassicalMultiDimensionalScaling,"CLASSICAL_MDS")
+
+void ClassicalMultiDimensionalScaling::registerKeywords( Keywords& keys ){
+ DimensionalityReductionBase::registerKeywords( keys );
+}
+
+ClassicalMultiDimensionalScaling::ClassicalMultiDimensionalScaling( const ActionOptions& ao):
+Action(ao),
+DimensionalityReductionBase(ao)
+{
+ if( dimredbase ) error("input to CLASSICAL_MDS should not be output from dimensionality reduction object");
+}
+
+void ClassicalMultiDimensionalScaling::calculateProjections( const Matrix<double>& targets, Matrix<double>& projections ){
+ // Retrieve the distances from the dimensionality reduction object
+ double half=(-0.5); Matrix<double> distances( half*targets );
+
+ // Apply centering transtion
+ unsigned n=distances.nrows(); double sum;
+ // First HM
+ for(unsigned i=0;i<n;++i){
+ sum=0; for(unsigned j=0;j<n;++j) sum+=distances(i,j);
+ for(unsigned j=0;j<n;++j) distances(i,j) -= sum/n;
+ }
+ // Now (HM)H
+ for(unsigned i=0;i<n;++i){
+ sum=0; for(unsigned j=0;j<n;++j) sum+=distances(j,i);
+ for(unsigned j=0;j<n;++j) distances(j,i) -= sum/n;
+ }
+
+ // Diagonalize matrix
+ std::vector<double> eigval(n); Matrix<double> eigvec(n,n);
+ diagMat( distances, eigval, eigvec );
+
+ // Pass final projections to map object
+ for(unsigned i=0;i<n;++i){
+ for(unsigned j=0;j<projections.ncols();++j) projections(i,j)=sqrt(eigval[n-1-j])*eigvec(n-1-j,i);
+ }
+}
+
+}
+}
diff --git a/src/dimred/DimensionalityReductionBase.cpp b/src/dimred/DimensionalityReductionBase.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f930cd75fbdfec83e3e7e722b42d0c976d8cfbd6
--- /dev/null
+++ b/src/dimred/DimensionalityReductionBase.cpp
@@ -0,0 +1,147 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2012 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.0.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "DimensionalityReductionBase.h"
+#include "reference/ReferenceConfiguration.h"
+#include "reference/MetricRegister.h"
+#include "core/PlumedMain.h"
+#include "core/Atoms.h"
+
+namespace PLMD {
+namespace dimred {
+
+void DimensionalityReductionBase::registerKeywords( Keywords& keys ){
+ analysis::AnalysisBase::registerKeywords( keys );
+ keys.add("compulsory","NLOW_DIM","number of low-dimensional coordinates required");
+}
+
+DimensionalityReductionBase::DimensionalityReductionBase( const ActionOptions& ao ):
+Action(ao),
+analysis::AnalysisBase(ao),
+myref(NULL),
+dimredbase(NULL)
+{
+ // Check that some dissimilarity information is available
+ if( !dissimilaritiesWereSet() ) error("dissimilarities have not been calcualted in input actions");
+ // Now we check if the input was a dimensionality reduction object
+ if( mydata ) dimredbase = dynamic_cast<DimensionalityReductionBase*>( mydata );
+
+ // Retrieve the dimension in the low dimensionality space
+ if( dimredbase ){
+ nlow=dimredbase->nlow;
+ } else {
+ parse("NLOW_DIM",nlow);
+ if( nlow<1 ) error("dimensionality of low dimensional space must be at least one");
+ }
+ log.printf(" projecting in %d dimensional space \n",nlow);
+
+ ReferenceConfigurationOptions("EUCLIDEAN");
+ myref=metricRegister().create<ReferenceConfiguration>("EUCLIDEAN");
+ std::vector<std::string> dimnames(nlow); std::string num;
+ for(unsigned i=0;i<nlow;++i){ Tools::convert(i+1,num); dimnames[i] = getLabel() + "." + num; }
+ myref->setNamesAndAtomNumbers( std::vector<AtomNumber>(), dimnames );
+}
+
+DimensionalityReductionBase::~DimensionalityReductionBase(){
+ delete myref;
+}
+
+ReferenceConfiguration* DimensionalityReductionBase::getReferenceConfiguration( const unsigned& idata ){
+ std::vector<double> pp(nlow); for(unsigned i=0;i<nlow;++i) pp[i]=projections(idata,i);
+ std::vector<double> empty( pp.size() );
+ myref->setReferenceConfig( std::vector<Vector>(), pp, empty );
+ return myref;
+}
+
+// double DimensionalityReductionBase::getInputDissimilarity( const unsigned& idata, const unsigned& jdata ){
+// if( dimredbase && use_dimred_dissims ) return dimredbase->getOutputDissimilarity( idata, jdata );
+// if( dimredbase ) return dimredbase->getInputDissimilarity( idata, jdata );
+// return getDissimilarity( idata, jdata );
+// }
+
+void DimensionalityReductionBase::getDataPoint( const unsigned& idata, std::vector<double>& point, double& weight ) const {
+ if( point.size()!=nlow ) point.resize( nlow );
+ weight = getWeight(idata); for(unsigned i=0;i<nlow;++i) point[i]=projections(idata,i);
+}
+
+// double DimensionalityReductionBase::getOutputDissimilarity( const unsigned& idata, const unsigned& jdata ){
+// double dissim=0; for(unsigned i=0;i<nlow;++i){ double tmp=projections(idata,i)-projections(jdata,i); dissim+=tmp*tmp; }
+// return dissim;
+// }
+
+void DimensionalityReductionBase::performAnalysis(){
+ // Resize the tempory array (this is used for out of sample)
+ dtargets.resize( getNumberOfDataPoints() );
+ // Resize the projections array
+ projections.resize( getNumberOfDataPoints(), nlow );
+
+// // Retrieve the weights from the previous calculation
+// std::vector<double> lweights( getNumberOfDataPoints() );
+// for(unsigned i=0;i<getNumberOfDataPoints();++i) lweights[i]=getWeight(i);
+// setOutputWeights( lweights );
+
+ // Retreive the projections from the previous calculation
+ if( dimredbase ){
+ std::vector<double> newp( nlow ); double w;
+ for(unsigned i=0;i<getNumberOfDataPoints();++i){
+ dimredbase->getDataPoint( i, newp, w ); plumed_dbg_assert( newp.size()==nlow );
+ for(unsigned j=0;j<nlow;++j) projections(i,j)=newp[j];
+ }
+ }
+ // Calculate matrix of dissimilarities
+ Matrix<double> targets( getNumberOfDataPoints(), getNumberOfDataPoints() );
+ for(unsigned i=1;i<getNumberOfDataPoints();++i){
+ for(unsigned j=0;j<i;++j) targets(i,j)=targets(j,i)=getDissimilarity( i, j );
+ }
+ // This calculates the projections of the points
+ calculateProjections( targets, projections );
+}
+
+double DimensionalityReductionBase::calculateStress( const std::vector<double>& p, std::vector<double>& d ){
+
+ // Zero derivative and stress accumulators
+ for(unsigned i=0;i<p.size();++i) d[i]=0.0;
+ double stress=0; std::vector<double> dtmp( p.size() );
+
+ // Now accumulate total stress on system
+ for(unsigned i=0;i<dtargets.size();++i){
+ if( dtargets[i]<epsilon ) continue ;
+
+ // Calculate distance in low dimensional space
+ double dd=0;
+ for(unsigned j=0;j<p.size();++j){ dtmp[j]=p[j]-projections(i,j); dd+=dtmp[j]*dtmp[j]; }
+ dd = sqrt(dd);
+
+ // Now do transformations and calculate differences
+ double ddiff = dd - dtargets[i];
+
+ // Calculate derivatives
+ double pref = 2.*getWeight(i) / dd;
+ for(unsigned j=0;j<p.size();++j) d[j] += pref*ddiff*dtmp[j];
+
+ // Accumulate the total stress
+ stress += getWeight(i)*ddiff*ddiff;
+ }
+ return stress;
+}
+
+}
+}
diff --git a/src/dimred/DimensionalityReductionBase.h b/src/dimred/DimensionalityReductionBase.h
new file mode 100644
index 0000000000000000000000000000000000000000..8543dfc8703798438f879333eec00af3f9d53f9c
--- /dev/null
+++ b/src/dimred/DimensionalityReductionBase.h
@@ -0,0 +1,78 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2012-2014 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_dimred_DimensionalityReductionBase_h
+#define __PLUMED_dimred_DimensionalityReductionBase_h
+
+#include "analysis/AnalysisBase.h"
+
+namespace PLMD {
+namespace dimred {
+
+class DimensionalityReductionBase : public analysis::AnalysisBase {
+friend class ProjectNonLandmarkPoints;
+friend class SketchMapBase;
+private:
+/// This are the target distances for a single point.
+/// This is used when we do out of sample or pointwise global optimization
+ std::vector<double> dtargets;
+/// We create a reference configuration here so that we can pass projection data
+/// quickly
+ ReferenceConfiguration* myref;
+/// The projections that were generated by the dimensionality reduction algorithm
+ Matrix<double> projections;
+protected:
+/// Dimensionality of low dimensional space
+ unsigned nlow;
+/// A pointer to any dimensionality reduction base that we have got projections from
+ DimensionalityReductionBase* dimredbase;
+public:
+ static void registerKeywords( Keywords& keys );
+ DimensionalityReductionBase( const ActionOptions& );
+ ~DimensionalityReductionBase();
+/// Get the ith data point (this returns the projection)
+ void getDataPoint( const unsigned& idata, std::vector<double>& point, double& weight ) const ;
+/// Get a reference configuration (this returns the projection)
+ ReferenceConfiguration* getReferenceConfiguration( const unsigned& idata );
+/// Actually perform the analysis
+ void performAnalysis();
+/// Calculate the projections of points
+ virtual void calculateProjections( const Matrix<double>& , Matrix<double>& )=0;
+/// Set one of the elements in the target vector. This target vector is used
+/// when we use calculateStress when finding the projections of individual points.
+/// For example this function is used in PLMD::dimred::ProjectOutOfSample
+ virtual void setTargetDistance( const unsigned& , const double& );
+/// Calculate the pointwise stress on one point when it is located at pp.
+/// This function makes use of the distance data in dtargets
+/// It is used in PLMD::dimred::ProjectOutOfSample and in pointwise optimisation
+ virtual double calculateStress( const std::vector<double>& pp, std::vector<double>& der );
+/// Overwrite virtual function in ActionWithVessel
+ void performTask( const unsigned& , const unsigned& , MultiValue& ) const { plumed_error(); }
+};
+
+inline
+void DimensionalityReductionBase::setTargetDistance( const unsigned& idata, const double& dist ){
+ dtargets[idata]=dist;
+}
+
+}
+}
+#endif
diff --git a/src/dimred/Makefile b/src/dimred/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..5d5eaab8eb4624680b2168b412eceef3da7c40eb
--- /dev/null
+++ b/src/dimred/Makefile
@@ -0,0 +1,4 @@
+USE=core tools reference analysis
+
+# generic makefile
+include ../maketools/make.module
diff --git a/src/dimred/ProjectNonLandmarkPoints.cpp b/src/dimred/ProjectNonLandmarkPoints.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5b039c47d053a3499b8ec8b9583fcb2e9ea72234
--- /dev/null
+++ b/src/dimred/ProjectNonLandmarkPoints.cpp
@@ -0,0 +1,129 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2012 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.0.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "core/ActionRegister.h"
+#include "core/PlumedMain.h"
+#include "core/ActionSet.h"
+#include "tools/Random.h"
+#include "reference/MetricRegister.h"
+#include "tools/ConjugateGradient.h"
+#include "analysis/AnalysisBase.h"
+#include "DimensionalityReductionBase.h"
+
+namespace PLMD {
+namespace dimred {
+
+class ProjectNonLandmarkPoints : public analysis::AnalysisBase {
+private:
+/// Tolerance for conjugate gradient algorithm
+ double cgtol;
+/// Number of diemsions in low dimensional space
+ unsigned nlow;
+/// We create a reference configuration here so that we can pass projection data
+/// quickly
+ ReferenceConfiguration* myref;
+/// The class that calcualtes the projection of the data that is required
+ DimensionalityReductionBase* mybase;
+/// Generate a projection of the ith data point - this is called in two routine
+ void generateProjection( const unsigned& idata, std::vector<double>& point );
+public:
+ static void registerKeywords( Keywords& keys );
+ ProjectNonLandmarkPoints( const ActionOptions& ao );
+ ~ProjectNonLandmarkPoints();
+/// Get the ith data point (this returns the projection)
+ void getDataPoint( const unsigned& idata, std::vector<double>& point );
+/// Get a reference configuration (this returns the projection)
+ ReferenceConfiguration* getReferenceConfiguration( const unsigned& idata );
+/// This does nothing -- projections are calculated when getDataPoint and getReferenceConfiguration are called
+ void performAnalysis(){}
+/// This just calls calculate stress in the underlying projection object
+ double calculateStress( const std::vector<double>& pp, std::vector<double>& der );
+/// Overwrite virtual function in ActionWithVessel
+ void performTask( const unsigned& , const unsigned& , MultiValue& ) const { plumed_error(); }
+};
+
+PLUMED_REGISTER_ACTION(ProjectNonLandmarkPoints,"PROJECT_ALL_ANALYSIS_DATA")
+
+void ProjectNonLandmarkPoints::registerKeywords( Keywords& keys ){
+ analysis::AnalysisBase::registerKeywords( keys );
+ keys.add("compulsory","PROJECTION","the projection that you wish to generate out-of-sample projections with");
+ keys.add("compulsory","CGTOL","1E-6","the tolerance for the conjugate gradient optimisation");
+}
+
+ProjectNonLandmarkPoints::ProjectNonLandmarkPoints( const ActionOptions& ao ):
+Action(ao),
+analysis::AnalysisBase(ao),
+mybase(NULL)
+{
+ std::string myproj; parse("PROJECTION",myproj);
+ mybase = plumed.getActionSet().selectWithLabel<DimensionalityReductionBase*>( myproj );
+ if( !mybase ) error("could not find projection of data named " + myproj );
+ nlow = mybase->nlow;
+
+ if( mybase->getBaseDataLabel()!=getBaseDataLabel() ) error("mismatch between base data labels");
+
+ log.printf(" generating out-of-sample projections using projection with label %s \n",myproj.c_str() );
+ parse("CGTOL",cgtol);
+
+ ReferenceConfigurationOptions("EUCLIDEAN");
+ myref=metricRegister().create<ReferenceConfiguration>("EUCLIDEAN");
+ std::vector<std::string> dimnames(nlow); std::string num;
+ for(unsigned i=0;i<nlow;++i){ Tools::convert(i+1,num); dimnames[i] = getLabel() + "." + num; }
+ myref->setNamesAndAtomNumbers( std::vector<AtomNumber>(), dimnames );
+}
+
+ProjectNonLandmarkPoints::~ProjectNonLandmarkPoints(){
+ delete myref;
+}
+
+void ProjectNonLandmarkPoints::generateProjection( const unsigned& idata, std::vector<double>& point ){
+ ConjugateGradient<ProjectNonLandmarkPoints> myminimiser( this );
+ unsigned closest=0; double mindist = sqrt( getDissimilarity( idata, mybase->getDataPointIndexInBase(0) ) );
+ mybase->setTargetDistance( 0, mindist );
+ for(unsigned i=1;i<mybase->getNumberOfDataPoints();++i){
+ double dist = sqrt( getDissimilarity( idata, mybase->getDataPointIndexInBase(i) ) );
+ mybase->setTargetDistance( i, dist );
+ if( dist<mindist ){ mindist=dist; closest=i; }
+ }
+ // Put the initial guess near to the closest landmark -- may wish to use grid here again Sandip??
+ Random random; random.setSeed(-1234);
+ for(unsigned j=0;j<nlow;++j) point[j]=mybase->projections(closest,j) + (random.RandU01() - 0.5)*0.01;
+ myminimiser.minimise( cgtol, point, &ProjectNonLandmarkPoints::calculateStress );
+}
+
+ReferenceConfiguration* ProjectNonLandmarkPoints::getReferenceConfiguration( const unsigned& idata ){
+ std::vector<double> pp(nlow); std::vector<double> empty( pp.size() ); generateProjection( idata, pp );
+ myref->setReferenceConfig( std::vector<Vector>(), pp, empty );
+ return myref;
+}
+
+void ProjectNonLandmarkPoints::getDataPoint( const unsigned& idata, std::vector<double>& point ){
+ if( point.size()!=nlow ) point.resize( nlow );
+ generateProjection( idata, point );
+}
+
+double ProjectNonLandmarkPoints::calculateStress( const std::vector<double>& pp, std::vector<double>& der ){
+ return mybase->calculateStress( pp, der );
+}
+
+}
+}
+
diff --git a/src/dimred/SMACOF.cpp b/src/dimred/SMACOF.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..aa748ce351a570ad6aefe89cf3bf4a829ae54cac
--- /dev/null
+++ b/src/dimred/SMACOF.cpp
@@ -0,0 +1,94 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2012 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.0.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "SMACOF.h"
+
+namespace PLMD {
+namespace dimred {
+
+void SMACOF::run( const Matrix<double>& Weights, const Matrix<double>& Distances, const double& tol, const unsigned& maxloops, Matrix<double>& InitialZ ){
+ unsigned M = Distances.nrows();
+
+ // Calculate V
+ Matrix<double> V(M,M); double totalWeight=0.;
+ for(unsigned i=0; i<M; ++i){
+ for(unsigned j=0; j<M; ++j){
+ if(i==j) continue;
+ V(i,j)=-Weights(i,j);
+ if( j<i ) totalWeight+=Weights(i,j);
+ }
+ for(unsigned j=0; j<M; ++j){
+ if(i==j)continue;
+ V(i,i)-=V(i,j);
+ }
+ }
+
+ // And pseudo invert V
+ Matrix<double> mypseudo(M, M); pseudoInvert(V, mypseudo);
+ Matrix<double> dists( M, M ); double myfirstsig = calculateSigma( Weights, Distances, InitialZ, dists ) / totalWeight;
+
+ // initial sigma is made up of the original distances minus the distances between the projections all squared.
+ Matrix<double> temp( M, M ), BZ( M, M ), newZ( M, InitialZ.ncols() );
+ for(unsigned n=0;n<maxloops;++n){
+ if(n==maxloops-1) plumed_merror("ran out of steps in SMACOF algorithm");
+
+ // Recompute BZ matrix
+ for(unsigned i=0; i<M; ++i){
+ for(unsigned j=0; j<M; ++j){
+ if(i==j) continue; //skips over the diagonal elements
+
+ if( dists(i,j)>0 ) BZ(i,j) = -Weights(i,j)*Distances(i,j) / dists(i,j);
+ else BZ(i,j)=0.;
+ }
+ //the diagonal elements are -off diagonal elements BZ(i,i)-=BZ(i,j) (Equation 8.25)
+ BZ(i,i)=0; //create the space memory for the diagonal elements which are scalars
+ for(unsigned j=0; j<M; ++j){
+ if(i==j) continue;
+ BZ(i,i)-=BZ(i,j);
+ }
+ }
+
+ mult( mypseudo, BZ, temp); mult(temp, InitialZ, newZ);
+ //Compute new sigma
+ double newsig = calculateSigma( Weights, Distances, newZ, dists ) / totalWeight;
+ printf("SIGMA VALUES %d %f %f %f \n",n, myfirstsig,newsig,myfirstsig-newsig);
+ //Computing whether the algorithm has converged (has the mass of the potato changed
+ //when we put it back in the oven!)
+ if( fabs( newsig - myfirstsig )<tol ) break;
+ myfirstsig=newsig;
+ InitialZ = newZ;
+ }
+}
+
+double SMACOF::calculateSigma( const Matrix<double>& Weights, const Matrix<double>& Distances, const Matrix<double>& InitialZ, Matrix<double>& dists ){
+ unsigned M = Distances.nrows(); double sigma=0;
+ for(unsigned i=1; i<M; ++i){
+ for(unsigned j=0; j<i; ++j){
+ double dlow=0; for(unsigned k=0; k<InitialZ.ncols();++k){ double tmp=InitialZ(i,k) - InitialZ(j,k); dlow+=tmp*tmp; }
+ dists(i,j)=dists(j,i)=sqrt(dlow); double tmp3 = Distances(i,j) - dists(i,j);
+ sigma += Weights(i,j)*tmp3*tmp3;
+ }
+ }
+ return sigma;
+}
+
+}
+}
diff --git a/src/dimred/SMACOF.h b/src/dimred/SMACOF.h
new file mode 100644
index 0000000000000000000000000000000000000000..4c79e91e6a65f078bdf290d1c53f0dca39793b1a
--- /dev/null
+++ b/src/dimred/SMACOF.h
@@ -0,0 +1,39 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2012 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.0.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_dimred_SMACOF_h
+#define __PLUMED_dimred_SMACOF_h
+
+#include <vector>
+#include "tools/Matrix.h"
+
+namespace PLMD {
+namespace dimred {
+
+class SMACOF {
+public:
+ static double calculateSigma( const Matrix<double>& Weights, const Matrix<double>& Distances, const Matrix<double>& InitialZ, Matrix<double>& dists );
+ static void run( const Matrix<double>& Weights, const Matrix<double>& Distances, const double& tol, const unsigned& maxloops, Matrix<double>& InitialZ);
+};
+
+}
+}
+#endif
diff --git a/src/dimred/SketchMapBase.cpp b/src/dimred/SketchMapBase.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..87c29c2370f61e023212d01085ac97e156a93ad7
--- /dev/null
+++ b/src/dimred/SketchMapBase.cpp
@@ -0,0 +1,128 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2012 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.0.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "SketchMapBase.h"
+
+namespace PLMD {
+namespace dimred {
+
+void SketchMapBase::registerKeywords( Keywords& keys ){
+ DimensionalityReductionBase::registerKeywords( keys );
+ keys.remove("NLOW_DIM");
+ keys.add("compulsory","HIGH_DIM_FUNCTION","as in input action","the parameters of the switching function in the high dimensional space");
+ keys.add("compulsory","LOW_DIM_FUNCTION","as in input action","the parameters of the switching function in the low dimensional space");
+ keys.add("compulsory","MIXPARAM","0.0","the ammount of the pure distances to mix into the stress function");
+}
+
+SketchMapBase::SketchMapBase( const ActionOptions& ao ):
+Action(ao),
+DimensionalityReductionBase(ao),
+smapbase(NULL)
+{
+ // Check if we have data from a input sketch-map object - we can reuse switching functions wahoo!!
+ smapbase = dynamic_cast<SketchMapBase*>( dimredbase );
+
+ // Read in the switching functions
+ std::string linput,hinput, errors;
+ parse("HIGH_DIM_FUNCTION",hinput);
+ if( hinput=="as in input action" ){
+ if( !smapbase ) error("high dimensional switching funciton has not been set - use HIGH_DIM_FUNCTION");
+ reuse_hd=true;
+ log.printf(" reusing high dimensional filter function defined in previous sketch-map action\n");
+ } else {
+ reuse_hd=false;
+ highdf.set(hinput,errors);
+ if(errors.length()>0) error(errors);
+ log.printf(" filter function for dissimilarities in high dimensional space has cutoff %s \n",highdf.description().c_str() );
+ }
+
+ parse("LOW_DIM_FUNCTION",linput);
+ if( linput=="as in input action" ){
+ if( !smapbase ) error("low dimensional switching funciton has not been set - use LOW_DIM_FUNCTION");
+ reuse_ld=true;
+ log.printf(" reusing low dimensional filter function defined in previous sketch-map action\n");
+ } else {
+ reuse_ld=false;
+ lowdf.set(hinput,errors);
+ if(errors.length()>0) error(errors);
+ log.printf(" filter function for distances in low dimensionality space has cutoff %s \n",highdf.description().c_str() );
+ }
+
+ // Read the mixing parameter
+ parse("MIXPARAM",mixparam);
+ if( mixparam<0 || mixparam>1 ) error("mixing parameter must be between 0 and 1");
+ log.printf(" mixing %f of pure distances with %f of filtered distances \n",mixparam,1.-mixparam);
+}
+
+void SketchMapBase::calculateProjections( const Matrix<double>& targets, Matrix<double>& projections ){
+ // These hold data so that we can do stress calculations
+ dtargets.resize( targets.nrows() ); ftargets.resize( targets.nrows() );
+
+ // Matrices for storing input data
+ Matrix<double> transformed( targets.nrows(), targets.ncols() );
+ Matrix<double> distances( targets.nrows(), targets.ncols() );
+
+ // Transform the high dimensional distances
+ double df; distances=0.; transformed=0.;
+ for(unsigned i=1;i<distances.ncols();++i){
+ for(unsigned j=0;j<i;++j){
+ distances(i,j)=distances(j,i)=sqrt( targets(i,j) );
+ transformed(i,j)=transformed(j,i)=transformHighDimensionalDistance( distances(i,j), df );
+ }
+ }
+ // And minimse
+ minimise( transformed, distances, projections );
+}
+
+double SketchMapBase::calculateStress( const std::vector<double>& p, std::vector<double>& d ){
+
+ // Zero derivative and stress accumulators
+ for(unsigned i=0;i<p.size();++i) d[i]=0.0;
+ double stress=0; std::vector<double> dtmp( p.size() );
+
+ // Now accumulate total stress on system
+ for(unsigned i=0;i<ftargets.size();++i){
+ if( dtargets[i]<epsilon ) continue ;
+
+ // Calculate distance in low dimensional space
+ double dd=0;
+ for(unsigned j=0;j<p.size();++j){ dtmp[j]=p[j]-projections(i,j); dd+=dtmp[j]*dtmp[j]; }
+ dd = sqrt(dd);
+
+ // Now do transformations and calculate differences
+ double df, fd = transformLowDimensionalDistance( dd, df );
+ double ddiff = dd - dtargets[i];
+ double fdiff = fd - ftargets[i];
+
+ // Calculate derivatives
+ double pref = 2.*getWeight(i) / dd ;
+ for(unsigned j=0;j<p.size();++j) d[j] += pref*( (1-mixparam)*fdiff*df + mixparam*ddiff )*dtmp[j];
+
+ // Accumulate the total stress
+ stress += getWeight(i)*( (1-mixparam)*fdiff*fdiff + mixparam*ddiff*ddiff );
+ }
+ return stress;
+}
+
+}
+}
+
+
diff --git a/src/dimred/SketchMapBase.h b/src/dimred/SketchMapBase.h
new file mode 100644
index 0000000000000000000000000000000000000000..1026beb211ebb6fbbcf147049de8988ccbecd53e
--- /dev/null
+++ b/src/dimred/SketchMapBase.h
@@ -0,0 +1,86 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2012-2014 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_dimred_SketchMapBase_h
+#define __PLUMED_dimred_SketchMapBase_h
+
+#include "tools/SwitchingFunction.h"
+#include "DimensionalityReductionBase.h"
+
+namespace PLMD {
+namespace dimred {
+
+class SketchMapBase : public DimensionalityReductionBase {
+private:
+/// To save us retyping switching functions many times the code will reuse the
+/// ones from previous sketch-map objects
+ bool reuse_hd, reuse_ld;
+/// Was previous action in chain a sketch-map action
+ SketchMapBase* smapbase;
+/// Switching functions for low and high dimensional space
+ SwitchingFunction lowdf, highdf;
+/// This is used within calculate stress to hold the target distances and the
+/// target values for the high dimensional switching function
+ std::vector<double> dtargets, ftargets;
+protected:
+/// The fraction of pure distances to mix in when optimising
+ double mixparam;
+public:
+ static void registerKeywords( Keywords& keys );
+ SketchMapBase( const ActionOptions& );
+/// This starts the process of calculating the projections
+ void calculateProjections( const Matrix<double>& , Matrix<double>& );
+/// This finishes the process of calculating the prjections
+ virtual void minimise( const Matrix<double>& , const Matrix<double>& , Matrix<double>& )=0;
+/// Apply the low dimensional switching function to the value val
+ double transformLowDimensionalDistance( const double& val, double& df ) const ;
+/// Apply the high dimensional switching function to the value val
+ double transformHighDimensionalDistance( const double& val, double& df ) const ;
+/// Set the target distance and from it calculate the target value for the switching function
+/// This target vector is used when we use calculateStress when finding the projections of individual points.
+/// For example this function is used in PLMD::dimred::ProjectOutOfSample
+ void setTargetDistance( const unsigned& idata, const double& dist );
+/// Calculate the pointwise stress on one point when it is located at p.
+/// This function makes use of the distance data in dtargets and ftargets
+/// It is used in PLMD::dimred::ProjectOutOfSample and in pointwise optimisation
+ double calculateStress( const std::vector<double>& p, std::vector<double>& d );
+};
+
+inline
+double SketchMapBase::transformLowDimensionalDistance( const double& val, double& df ) const {
+ if( reuse_ld ) return smapbase->transformLowDimensionalDistance( val, df );
+ double ans=1.0 - lowdf.calculate( val, df ); df*=-val; return ans;
+}
+
+inline
+double SketchMapBase::transformHighDimensionalDistance( const double& val, double& df ) const {
+ if( reuse_hd ) return smapbase->transformHighDimensionalDistance( val, df );
+ double ans=1.0 - highdf.calculate( val, df ); df*=-val; return ans;
+}
+
+inline
+void SketchMapBase::setTargetDistance( const unsigned& idata, const double& dist ){
+ double df; dtargets[idata]=dist; ftargets[idata]=transformHighDimensionalDistance( dist, df );
+}
+
+}
+}
+#endif
diff --git a/src/dimred/SketchMapPointwise.cpp b/src/dimred/SketchMapPointwise.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b9f18fca5b31245af33743504a2469cf643a02cb
--- /dev/null
+++ b/src/dimred/SketchMapPointwise.cpp
@@ -0,0 +1,112 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2012 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.0.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "core/ActionRegister.h"
+#include "SketchMapBase.h"
+#include "tools/ConjugateGradient.h"
+#include "tools/GridSearch.h"
+
+
+namespace PLMD {
+namespace dimred {
+
+class SketchMapPointwise : public SketchMapBase {
+private:
+ unsigned ncycles;
+ bool nointerpolation;
+ double cgtol, gbuf;
+ std::vector<unsigned> npoints;
+public:
+ static void registerKeywords( Keywords& keys );
+ SketchMapPointwise( const ActionOptions& ao );
+ void minimise( const Matrix<double>& , const Matrix<double>& , Matrix<double>& );
+};
+
+PLUMED_REGISTER_ACTION(SketchMapPointwise,"SKETCHMAP_POINTWISE")
+
+void SketchMapPointwise::registerKeywords( Keywords& keys ){
+ SketchMapBase::registerKeywords( keys );
+ keys.add("compulsory","NCYCLES","5","the number of cycles of global optimisation to attempt");
+ keys.add("compulsory","CGTOL","1E-6","the tolerance for the conjugate gradient minimisation");
+ keys.add("compulsory","BUFFER","1.1","grid extent for search is (max projection - minimum projection) multiplied by this value");
+ keys.add("compulsory","NGRID","10","number of points to use in each grid direction");
+ keys.addFlag("NOINTERPOLATION",false,"do not use any interpolation to make the optimisation faster");
+}
+
+SketchMapPointwise::SketchMapPointwise( const ActionOptions& ao ):
+Action(ao),
+SketchMapBase(ao),
+npoints(nlow)
+{
+ parseVector("NGRID",npoints); parse("BUFFER",gbuf);
+ if( npoints.size()!=nlow ) error("vector giving number of grid point in each direction has wrong size");
+ parse("NCYCLES",ncycles); parse("CGTOL",cgtol);
+ parseFlag("NOINTERPOLATION",nointerpolation);
+
+ log.printf(" doing %d cycles of global optimisation sweeps\n",ncycles);
+ log.printf(" using grid of points that is %d",npoints[0]);
+ for(unsigned j=1;j<npoints.size();++j) log.printf(" by %d",npoints[j]);
+ log.printf(" and that is %f larger than the difference between the position of the minimum and maximum projection \n",gbuf);
+ log.printf(" tolerance for conjugate gradient algorithm equals %f \n",cgtol);
+ if( nointerpolation ) log.printf(" not using any interpolation in grid search\n");
+}
+
+void SketchMapPointwise::minimise( const Matrix<double>& transformed, const Matrix<double>& distances, Matrix<double>& projections ){
+ std::vector<double> gmin( nlow ), gmax( nlow ), mypoint( nlow );
+
+ // Find the extent of the grid
+ for(unsigned j=0;j<nlow;++j) gmin[j]=gmax[j]=projections(0,j);
+ for(unsigned i=1;i<getNumberOfDataPoints();++i){
+ for(unsigned j=0;j<nlow;++j){
+ if( projections(i,j) < gmin[j] ) gmin[j] = projections(i,j);
+ if( projections(i,j) > gmax[j] ) gmax[j] = projections(i,j);
+ }
+ }
+ for(unsigned j=0;j<nlow;++j){
+ double gbuffer = 0.5*gbuf*( gmax[j]-gmin[j] ) - 0.5*( gmax[j]- gmin[j] );
+ gmin[j]-=gbuffer; gmax[j]+=gbuffer;
+ }
+
+ // And do the search
+ ConjugateGradient<SketchMapPointwise> mycgminimise( this );
+ GridSearch<SketchMapPointwise> mygridsearch( gmin, gmax, npoints, this );
+ // Run multiple loops over all projections
+ for(unsigned i=0;i<ncycles;++i){
+ for(unsigned j=0;j<getNumberOfDataPoints();++j){
+ // Setup target distances and target functions for calculate stress
+ for(unsigned k=0;k<getNumberOfDataPoints();++k) setTargetDistance( k, distances(j,k) );
+
+ if( nointerpolation ){
+ // Minimise using grid search
+ mygridsearch.minimise( mypoint, &SketchMapPointwise::calculateStress );
+ // Minimise output using conjugate gradient
+ mycgminimise.minimise( cgtol, mypoint, &SketchMapPointwise::calculateStress );
+ } else {
+ plumed_merror("use a class here that Sandip will implement that allows you to do this with some interpolation");
+ }
+ // Save new projection
+ for(unsigned k=0;k<nlow;++k) projections(j,k)=mypoint[k];
+ }
+ }
+}
+
+}
+}
diff --git a/src/dimred/SketchMapSmacof.cpp b/src/dimred/SketchMapSmacof.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..591df4e5ce5eb10a23facba3e1710fd854d7616c
--- /dev/null
+++ b/src/dimred/SketchMapSmacof.cpp
@@ -0,0 +1,104 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2012 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.0.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "core/ActionRegister.h"
+#include "SketchMapBase.h"
+#include "SMACOF.h"
+
+namespace PLMD {
+namespace dimred {
+
+class SketchMapSmacof : public SketchMapBase {
+private:
+ unsigned max_smap, maxiter;
+ double smap_tol, iter_tol, regulariser;
+ double recalculateWeights( const Matrix<double>& transformed, const Matrix<double>& distances,
+ const Matrix<double>& projections, Matrix<double>& weights );
+public:
+ static void registerKeywords( Keywords& keys );
+ SketchMapSmacof( const ActionOptions& ao );
+ void minimise( const Matrix<double>& , const Matrix<double>& , Matrix<double>& );
+};
+
+PLUMED_REGISTER_ACTION(SketchMapSmacof,"SKETCHMAP_SMACOF")
+
+void SketchMapSmacof::registerKeywords( Keywords& keys ){
+ SketchMapBase::registerKeywords( keys );
+ keys.add("compulsory","SMACOF_TOL","1E-4","the tolerance for each SMACOF cycle");
+ keys.add("compulsory","SMACOF_MAXCYC","1000","maximum number of optimization cycles for SMACOF algorithm");
+ keys.add("compulsory","SMAP_TOL","1E-4","the tolerance for sketch-map");
+ keys.add("compulsory","SMAP_MAXCYC","100","maximum number of optimization cycles for iterative sketch-map algorithm");
+ keys.add("compulsory","REGULARISE_PARAM","0.001","this is used to ensure that we don't divide by zero when updating weights");
+}
+
+SketchMapSmacof::SketchMapSmacof( const ActionOptions& ao ):
+Action(ao),
+SketchMapBase(ao)
+{
+ parse("REGULARISE_PARAM",regulariser);
+ parse("SMACOF_MAXCYC",max_smap); parse("SMAP_MAXCYC",maxiter);
+ parse("SMACOF_TOL",smap_tol); parse("SMAP_TOL",iter_tol);
+}
+
+void SketchMapSmacof::minimise( const Matrix<double>& transformed, const Matrix<double>& distances, Matrix<double>& projections ){
+ Matrix<double> weights( distances.nrows(), distances.ncols() ); weights=0.;
+ double filt = recalculateWeights( transformed, distances, projections, weights );
+
+ for(unsigned i=0;i<maxiter;++i){
+ SMACOF::run( weights, distances, smap_tol, max_smap, projections );
+ // Recalculate weights matrix and sigma
+ double newsig = recalculateWeights( transformed, distances, projections, weights );
+ printf("HELLO SMACOF %d %f %f \n",i,newsig,filt);
+ // Test whether or not the algorithm has converged
+ if( fabs( newsig - filt )<iter_tol ) break;
+ // Make initial sigma into new sigma so that the value of new sigma is used every time so that the error can be reduced
+ filt=newsig;
+ }
+}
+
+double SketchMapSmacof::recalculateWeights( const Matrix<double>& transformed, const Matrix<double>& distances,
+ const Matrix<double>& projections, Matrix<double>& weights ){
+ double filt=0, totalWeight=0.;; double dr;
+ for(unsigned i=1; i<weights.nrows(); ++i){
+ for(unsigned j=0; j<i; ++j){
+ double ninj=getWeight(i)*getWeight(j); totalWeight += ninj;
+
+ double tempd=0;
+ for(unsigned k=0;k<projections.ncols();++k){
+ double tmp = projections(i,k) - projections(j,k);
+ tempd += tmp*tmp;
+ }
+ double dij=sqrt(tempd);
+
+ double fij = transformLowDimensionalDistance( dij, dr );
+ double filter=transformed(i,j)-fij;
+ double diff=distances(i,j) - dij;
+
+ if( fabs(diff)<regulariser ) weights(i,j)=weights(j,i)=0.0;
+ else weights(i,j)=weights(j,i) = ninj*( (1-mixparam)*( filter*dr )/diff + mixparam );
+ filt += ninj*( (1-mixparam)*filter*filter + mixparam*diff*diff );
+ }
+ }
+ return filt / totalWeight;
+}
+
+}
+}
diff --git a/src/dimred/SmacoffMDS.cpp b/src/dimred/SmacoffMDS.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3c41163ebb52adac9cab0cf7f0b227e2e7c56946
--- /dev/null
+++ b/src/dimred/SmacoffMDS.cpp
@@ -0,0 +1,84 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2012-2014 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "DimensionalityReductionBase.h"
+#include "core/ActionRegister.h"
+#include "SMACOF.h"
+
+//+PLUMEDOC ANALYSIS SMACOFF_MDS
+/*
+Create a low-dimensional projection of a trajectory using the classical multidimensional
+scaling algorithm.
+
+\par Examples
+
+*/
+//+ENDPLUMEDOC
+
+namespace PLMD {
+namespace dimred {
+
+class SmacofMDS : public DimensionalityReductionBase {
+private:
+ unsigned maxloops;
+ double tol;
+public:
+ static void registerKeywords( Keywords& keys );
+ SmacofMDS( const ActionOptions& );
+ void calculateProjections( const Matrix<double>& , Matrix<double>& );
+};
+
+PLUMED_REGISTER_ACTION(SmacofMDS,"SMACOF_MDS")
+
+void SmacofMDS::registerKeywords( Keywords& keys ){
+ DimensionalityReductionBase::registerKeywords( keys );
+ keys.remove("NLOW_DIM");
+ keys.add("compulsory","SMACOF_TOL","1E-4","tolerance for the SMACOF optimization algorith");
+ keys.add("compulsory","SMACOF_MAXCYC","1000","maximum number of optimization cycles for SMACOF algorithm");
+}
+
+SmacofMDS::SmacofMDS( const ActionOptions& ao):
+Action(ao),
+DimensionalityReductionBase(ao)
+{
+ if( !dimredbase ) error("SMACOF must be initialised using output from dimensionality reduction object");
+
+ parse("SMACOF_TOL",tol); parse("SMACOF_MAXCYC",maxloops);
+ log.printf(" running smacof to convergence at %f or for a maximum of %d steps \n",tol,maxloops);
+}
+
+void SmacofMDS::calculateProjections( const Matrix<double>& targets, Matrix<double>& projections ){
+abort();
+ // Take the square root of all the distances and the weights
+ Matrix<double> weights( targets.nrows(), targets.ncols() );
+ Matrix<double> distances( targets.nrows(), targets.ncols() );
+ for(unsigned i=1;i<distances.ncols();++i){
+ for(unsigned j=0;j<i;++j){
+ distances(i,j)=distances(j,i)=sqrt( targets(i,j) );
+ weights(i,j)=weights(j,i)=getWeight(i)*getWeight(j);
+ }
+ }
+ // And run SMACOF
+ SMACOF::run( weights, targets, tol, maxloops, projections );
+}
+
+}
+}
diff --git a/src/dimred/module.type b/src/dimred/module.type
new file mode 100644
index 0000000000000000000000000000000000000000..de832730330473a1870bd368868640da677460ae
--- /dev/null
+++ b/src/dimred/module.type
@@ -0,0 +1 @@
+default-off
diff --git a/src/tools/ConjugateGradient.h b/src/tools/ConjugateGradient.h
new file mode 100644
index 0000000000000000000000000000000000000000..68281a3dd60b176e6c8030097feefb72a01b37ea
--- /dev/null
+++ b/src/tools/ConjugateGradient.h
@@ -0,0 +1,66 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2013 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_tools_ConjugateGradient_h
+#define __PLUMED_tools_ConjugateGradient_h
+
+#include "MinimiseBase.h"
+
+namespace PLMD{
+
+template <class FCLASS>
+class ConjugateGradient : public MinimiseBase<FCLASS> {
+private:
+/// This is the pointer to the member funciton in the energy
+/// calculating class that calculates the energy
+ typedef double(FCLASS::*engf_pointer)( const std::vector<double>& p, std::vector<double>& der );
+ const unsigned ITMAX;
+ const double EPS;
+public:
+ ConjugateGradient( FCLASS* funcc ) : MinimiseBase<FCLASS>(funcc), ITMAX(200), EPS(1E-10) {}
+ void minimise( const double& ftol, std::vector<double>& p, engf_pointer myfunc );
+};
+
+template <class FCLASS>
+void ConjugateGradient<FCLASS>::minimise( const double& ftol, std::vector<double>& p, engf_pointer myfunc ){
+ std::vector<double> xi( p.size() ), g( p.size() ), h( p.size() );
+ double fp = this->calcDerivatives( p, xi, myfunc );
+ for(unsigned j=0;j<p.size();++j){ g[j] = -xi[j]; xi[j]=h[j]=g[j]; }
+
+ for(unsigned its=0;its<ITMAX;++its){
+ double fret=this->linemin( xi, p, myfunc );
+ // The exit condition
+ if( 2.0*fabs(fret-fp) <= ftol*(fabs(fret)+fabs(fp)+EPS)) { return; }
+ fp = fret; double igeng = this->calcDerivatives( p, xi, myfunc );
+ double ddg=0., gg=0.;
+ for(unsigned j=0;j<p.size();++j){ gg += g[j]*g[j]; ddg += (xi[j]+g[j])*xi[j]; }
+
+ if( gg==0.0 ) return;
+
+ double gam=ddg/gg;
+ for(unsigned j=0;j<p.size();++j){ g[j] = -xi[j]; xi[j]=h[j]=g[j]+gam*h[j]; }
+ }
+ plumed_merror("Too many interactions in conjugate gradient");
+}
+
+}
+
+#endif
diff --git a/src/tools/GridSearch.h b/src/tools/GridSearch.h
new file mode 100644
index 0000000000000000000000000000000000000000..7d52ab226c8671b8cc38a5aef80707a6991a699b
--- /dev/null
+++ b/src/tools/GridSearch.h
@@ -0,0 +1,85 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2013 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_tools_GridSearch_h
+#define __PLUMED_tools_GridSearch_h
+
+#include "MinimiseBase.h"
+#include <iostream>
+#include <math.h>
+namespace PLMD{
+
+template <class FCLASS>
+class GridSearch {
+private:
+/// This is the pointer to the member funciton in the energy
+/// calculating class that calculates the energy
+ typedef double(FCLASS::*engf_pointer)( const std::vector<double>& p, std::vector<double>& der );
+ FCLASS* myclass_func;
+ std::vector<unsigned> ngrid;
+ std::vector<double> min, delr;
+ unsigned np;
+public:
+ GridSearch( const std::vector<double>& mmin, const std::vector<double>& mmax, const std::vector<unsigned>& ng, FCLASS* funcc ) :
+ myclass_func( funcc ),
+ ngrid(ng),
+ min(mmin),
+ delr(mmin.size())
+ {
+ // Work out the stride
+ for(unsigned i=0;i<delr.size();++i) delr[i] = ( mmax[i] - mmin[i] ) / static_cast<double>( ng[i] );
+ // And the number of poitns in the grid
+ np=1; for(unsigned i=0;i<ngrid.size();++i) np*=ngrid[i];
+ }
+ void minimise( std::vector<double>& p, engf_pointer myfunc );
+};
+
+template <class FCLASS>
+void GridSearch<FCLASS>::minimise( std::vector<double>& p, engf_pointer myfunc ){
+
+ std::vector<double> fake_der( p.size() ); std::vector<unsigned> ind( p.size() );
+ unsigned pmin=0; double emin=(myclass_func->*myfunc)( min, fake_der );
+ for(unsigned i=1;i<np;++i){
+ unsigned myind=i; ind[0]=myind%ngrid[0];
+ for(unsigned j=1;j<p.size()-1;++j){
+ myind=(myind-ind[j-1])/ngrid[j-1];
+ ind[j]=myind%ngrid[j];
+ }
+ if( p.size()>=2 ) ind[p.size()-1] = (myind - ind[p.size()-2])/ngrid[p.size()-2];
+ for(unsigned j=0;j<p.size();++j) p[j] = min[j] + ind[j]*delr[j];
+
+ double eng = (myclass_func->*myfunc)( p, fake_der );
+ if( eng<emin ){ emin=eng; pmin=i; }
+ }
+
+ // Now recover lowest energy point
+ ind[0]=pmin%ngrid[0];
+ for(unsigned j=1;j<p.size()-1;++j){
+ pmin=(pmin-ind[j-1])/ngrid[j-1];
+ ind[j]=pmin%ngrid[j];
+ }
+ if( p.size()>=2 ) ind[p.size()-1] = (pmin - ind[p.size()-2])/ngrid[p.size()-2];
+ for(unsigned j=0;j<p.size();++j) p[j] = min[j] + ind[j]*delr[j];
+}
+
+}
+#endif
+
diff --git a/src/tools/Minimise1DBrent.h b/src/tools/Minimise1DBrent.h
new file mode 100644
index 0000000000000000000000000000000000000000..7339a53195c50e9a11fad7991d092d216db9d7e2
--- /dev/null
+++ b/src/tools/Minimise1DBrent.h
@@ -0,0 +1,189 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2013 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_tools_Minimise1DBrent_h
+#define __PLUMED_tools_Minimise1DBrent_h
+
+#include <vector>
+#include <string>
+
+namespace PLMD{
+
+/// A class for doing parabolic interpolation and minimisation of
+/// 1D functions using Brent's method.
+template <class FCLASS>
+class Minimise1DBrent {
+private:
+/// Has the minimum been bracketed
+ bool bracketed;
+/// Has the function been minimised
+ bool minimised;
+/// The tolerance for the line minimiser
+ double tol;
+/// The default ration by which successive intervals are magnified
+ const double GOLD;
+/// The maximum magnification allowed for a parabolic fit step
+ const double GLIMIT;
+/// Use to prevent any possible division by zero
+ const double TINY;
+/// Maximum number of interactions in line minimiser
+ const int ITMAX;
+/// The value of the golden ratio
+ const double CGOLD;
+/// A small number that protects against trying to achieve fractional
+/// accuracy for a minimum that happens to be exactly zero
+ const double ZEPS;
+/// This is the type specifier for the function to minimise
+ typedef double(FCLASS::*eng_pointer)( const double& val );
+/// Three points bracketting the minimum and the corresponding function values
+ double ax,bx,cx,fa,fb,fc,fmin;
+/// The class containing the function we are trying to minimise
+ FCLASS myclass_func;
+public:
+ Minimise1DBrent( const FCLASS& pf, const double& t=3.0E-8 );
+/// Bracket the minium
+ void bracket( const double& ax, const double& xx, eng_pointer eng );
+/// Find the minimum between two brackets
+ double minimise( eng_pointer eng );
+/// Return the value of the function at the minimum
+ double getMinimumValue() const ;
+};
+
+template <class FCLASS>
+Minimise1DBrent<FCLASS>::Minimise1DBrent( const FCLASS& pf, const double& t ):
+bracketed(false),
+tol(t),
+GOLD(1.618034),
+GLIMIT(100.0),
+TINY(1.0E-20),
+ITMAX(100),
+CGOLD(0.3819660),
+ZEPS(epsilon*1.0E-3),
+myclass_func(pf)
+{
+}
+
+template <class FCLASS>
+void Minimise1DBrent<FCLASS>::bracket( const double& a, const double& b, eng_pointer eng ){
+ ax=a; bx=b; double fu;
+ fa=(myclass_func.*eng)(ax); fb=(myclass_func.*eng)(bx);
+ if( fb>fa ){
+ double tmp;
+ tmp=ax; ax=bx; bx=tmp;
+ tmp=fa; fa=fb; fb=tmp;
+ }
+ cx=bx+GOLD*(bx-ax);
+ fc=(myclass_func.*eng)(cx);
+ while ( fb > fc ){
+ double r=(bx-ax)*(fb-fc);
+ double q=(bx-cx)*(fb-fa);
+ double u=bx-((bx-cx)*q-(bx-ax)*r)/(2.0*(fabs(q-r)>TINY?fabs(q-r):TINY)*(q-r>=0?1:-1));
+ double ulim=bx+GLIMIT*(cx-bx);
+ if((bx-u)*(u-cx) > 0.0 ){
+ fu=(myclass_func.*eng)(u);
+ if( fu < fc ){
+ ax=bx; bx=u; fa=fb; fb=fu; bracketed=true; return;
+ } else if( fu > fb ){
+ cx=u; fc=fu; bracketed=true; return;
+ }
+ u=cx+GOLD*(cx-bx); fu=(myclass_func.*eng)(u);
+ } else if((cx-u)*(u-ulim) > 0.0 ){
+ fu=(myclass_func.*eng)(u);
+ if( fu<fc ){
+ bx=cx; cx=u; u+=GOLD*(u-bx);
+ fb=fc; fc=fu; fu=(myclass_func.*eng)(u);
+ }
+ } else if( (u-ulim)*(ulim-cx) >= 0.0 ){
+ u=ulim;
+ fu=(myclass_func.*eng)(u);
+ } else {
+ u=cx+GOLD*(cx-bx);
+ fu=(myclass_func.*eng)(u);
+ }
+ ax=bx; bx=cx; cx=u;
+ fa=fb; fb=fc; fc=fu;
+ }
+ bracketed=true;
+}
+
+template <class FCLASS>
+double Minimise1DBrent<FCLASS>::minimise( eng_pointer eng ){
+ plumed_dbg_assert( bracketed );
+
+ double a,b,d=0.0,etemp,fu,fv,fw,fx;
+ double p,q,r,tol1,tol2,u,v,w,x,xm;
+ double e=0.0;
+
+ a=(ax < cx ? ax : cx );
+ b=(ax >= cx ? ax : cx );
+ x=w=v=bx;
+ fw=fv=fx=(myclass_func.*eng)(x);
+ for(unsigned iter=0;iter<ITMAX;++iter){
+ xm=0.5*(a+b);
+ tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
+ if( fabs(x-xm) <= (tol2-0.5*(b-a))) {
+ fmin=fx; minimised=true; return x;
+ }
+ if( fabs(e) > tol1 ){
+ r=(x-w)*(fx-fv);
+ q=(x-v)*(fx-fw);
+ p=(x-v)*q-(x-w)*r;
+ q=2.0*(q-r);
+ if( q > 0.0 ) p = -p;
+ q=fabs(q);
+ etemp=e;
+ e=d;
+ if( fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x) ){
+ d = CGOLD*(e=(x >= xm ? a-x : b-x ));
+ } else {
+ d=p/q; u=x+d;
+ if(u-a < tol2 || b-u < tol2 ) d=(xm-x>=0?fabs(tol1):-fabs(tol1));
+ }
+ } else {
+ d=CGOLD*(e=( x >= xm ? a-x : b-x ));
+ }
+ if( fabs(d)>=tol1) u=x+d; else u=x+(d>=0?fabs(tol1):-fabs(tol1));
+ fu=(myclass_func.*eng)(u);
+ if( fu <= fx ){
+ if( u >= x ) a=x; else b=x;
+ v=w; fv=fw;
+ w=x; fw=fx;
+ x=u; fx=fu;
+ } else {
+ if( u < x ) a=u; else b=u;
+ if( fu <=fw || w==x ){
+ v=w; w=u; fv=fw; fw=fu;
+ } else if( fu <= fv || v==x || v==w ){
+ v=u; fv=fu;
+ }
+ }
+ }
+ plumed_merror("Too many interactions in brent");
+}
+
+template <class FCLASS>
+double Minimise1DBrent<FCLASS>::getMinimumValue() const {
+ plumed_dbg_assert( minimised );
+ return fmin;
+}
+
+}
+#endif
diff --git a/src/tools/MinimiseBase.h b/src/tools/MinimiseBase.h
new file mode 100644
index 0000000000000000000000000000000000000000..a6c62257352663a8fe0b7353d1c4c6e17cd5274d
--- /dev/null
+++ b/src/tools/MinimiseBase.h
@@ -0,0 +1,111 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2013 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_tools_MinimiseBase_h
+#define __PLUMED_tools_MinimiseBase_h
+
+#include "Minimise1DBrent.h"
+
+namespace PLMD {
+
+template <class FCLASS>
+class F1dim {
+private:
+/// This is the pointer to the member funciton in the energy
+/// calculating class that calculates the energy
+ typedef double(FCLASS::*engf_pointer)( const std::vector<double>& p, std::vector<double>& der );
+/// Pointer to the vector containing an initial position on the vector
+ const std::vector<double>& p;
+/// The direction of the vector we are minimising along
+ const std::vector<double>& dir;
+/// Tempory vector that holds a point at which we want to calculate the energy
+ std::vector<double> pt;
+/// Vector that holds the derivatives at the point at which we calculate the energy (these are not used)
+ std::vector<double> fake_der;
+/// Class containging the function in the class
+ FCLASS* func;
+/// Member of class that calculates the energy we are trying to mnimise
+ engf_pointer calc;
+public:
+ F1dim( const std::vector<double>& pp, const std::vector<double>& dd, FCLASS* ff, engf_pointer cc );
+/// Calculate the energy at \f$\mathbf{p} + xt*\mathbf{dir}\f$
+ double getEng( const double& xt );
+};
+
+template <class FCLASS>
+F1dim<FCLASS>::F1dim( const std::vector<double>& pp, const std::vector<double>& dd, FCLASS* ff, engf_pointer cc ):
+p(pp),
+dir(dd),
+pt(pp.size()),
+fake_der(pp.size()),
+func(ff),
+calc(cc)
+{
+}
+
+template <class FCLASS>
+double F1dim<FCLASS>::getEng( const double& xt ){
+ for(unsigned j=0;j<pt.size();++j) pt[j] = p[j] + xt*dir[j];
+ return (func->*calc)(pt,fake_der);
+}
+
+template <class FCLASS>
+class MinimiseBase {
+private:
+/// This is the pointer to the member funciton in the energy
+/// calculating class that calculates the energy
+ typedef double(FCLASS::*engf_pointer)( const std::vector<double>& p, std::vector<double>& der );
+/// The class that calculates the energy given a position
+ FCLASS* myclass_func;
+/// Member of class that calculates the energy we are trying to minimise
+ engf_pointer calc;
+protected:
+/// This is the line minimiser
+ double linemin( const std::vector<double>& dir, std::vector<double>& p, engf_pointer myfunc );
+/// This calculates the derivatives at a point
+ double calcDerivatives( const std::vector<double>& p, std::vector<double>& der, engf_pointer myfunc );
+public:
+ MinimiseBase( FCLASS* funcc ) : myclass_func(funcc) {}
+};
+
+template <class FCLASS>
+double MinimiseBase<FCLASS>::linemin( const std::vector<double>& dir, std::vector<double>& p, engf_pointer myfunc ){
+ // Construct the object that turns points on a line into vectors
+ F1dim<FCLASS> f1dim( p, dir, myclass_func, myfunc );
+
+ // Construct an object that will do the line search for the minimum
+ Minimise1DBrent<F1dim<FCLASS> > bb(f1dim);
+
+ // This does the actual line minimisation
+ double ax=0.0, xx=1.0;
+ bb.bracket( ax, xx, &F1dim<FCLASS>::getEng );
+ double xmin=bb.minimise( &F1dim<FCLASS>::getEng );
+ for(unsigned i=0;i<p.size();++i) p[i] += xmin*dir[i];
+ return bb.getMinimumValue();
+}
+
+template <class FCLASS>
+double MinimiseBase<FCLASS>::calcDerivatives( const std::vector<double>& p, std::vector<double>& der, engf_pointer myfunc ){
+ return (myclass_func->*myfunc)( p, der );
+}
+
+}
+#endif