1 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 Copyright (c) 2015-2021 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 "SMACOF.h"
23
24 namespace PLMD {
25 namespace dimred {
26
run(const Matrix<double> & Weights,const Matrix<double> & Distances,const double & tol,const unsigned & maxloops,Matrix<double> & InitialZ)27 void SMACOF::run( const Matrix<double>& Weights, const Matrix<double>& Distances, const double& tol, const unsigned& maxloops, Matrix<double>& InitialZ ) {
28 unsigned M = Distances.nrows();
29
30 // Calculate V
31 Matrix<double> V(M,M); double totalWeight=0.;
32 for(unsigned i=0; i<M; ++i) {
33 for(unsigned j=0; j<M; ++j) {
34 if(i==j) continue;
35 V(i,j)=-Weights(i,j);
36 if( j<i ) totalWeight+=Weights(i,j);
37 }
38 for(unsigned j=0; j<M; ++j) {
39 if(i==j)continue;
40 V(i,i)-=V(i,j);
41 }
42 }
43
44 // And pseudo invert V
45 Matrix<double> mypseudo(M, M); pseudoInvert(V, mypseudo);
46 Matrix<double> dists( M, M ); double myfirstsig = calculateSigma( Weights, Distances, InitialZ, dists ) / totalWeight;
47
48 // initial sigma is made up of the original distances minus the distances between the projections all squared.
49 Matrix<double> temp( M, M ), BZ( M, M ), newZ( M, InitialZ.ncols() );
50 for(unsigned n=0; n<maxloops; ++n) {
51 if(n==maxloops-1) plumed_merror("ran out of steps in SMACOF algorithm");
52
53 // Recompute BZ matrix
54 for(unsigned i=0; i<M; ++i) {
55 for(unsigned j=0; j<M; ++j) {
56 if(i==j) continue; //skips over the diagonal elements
57
58 if( dists(i,j)>0 ) BZ(i,j) = -Weights(i,j)*Distances(i,j) / dists(i,j);
59 else BZ(i,j)=0.;
60 }
61 //the diagonal elements are -off diagonal elements BZ(i,i)-=BZ(i,j) (Equation 8.25)
62 BZ(i,i)=0; //create the space memory for the diagonal elements which are scalars
63 for(unsigned j=0; j<M; ++j) {
64 if(i==j) continue;
65 BZ(i,i)-=BZ(i,j);
66 }
67 }
68
69 mult( mypseudo, BZ, temp); mult(temp, InitialZ, newZ);
70 //Compute new sigma
71 double newsig = calculateSigma( Weights, Distances, newZ, dists ) / totalWeight;
72 //Computing whether the algorithm has converged (has the mass of the potato changed
73 //when we put it back in the oven!)
74 if( fabs( newsig - myfirstsig )<tol ) break;
75 myfirstsig=newsig;
76 InitialZ = newZ;
77 }
78 }
79
calculateSigma(const Matrix<double> & Weights,const Matrix<double> & Distances,const Matrix<double> & InitialZ,Matrix<double> & dists)80 double SMACOF::calculateSigma( const Matrix<double>& Weights, const Matrix<double>& Distances, const Matrix<double>& InitialZ, Matrix<double>& dists ) {
81 unsigned M = Distances.nrows(); double sigma=0;
82 for(unsigned i=1; i<M; ++i) {
83 for(unsigned j=0; j<i; ++j) {
84 double dlow=0; for(unsigned k=0; k<InitialZ.ncols(); ++k) { double tmp=InitialZ(i,k) - InitialZ(j,k); dlow+=tmp*tmp; }
85 dists(i,j)=dists(j,i)=sqrt(dlow); double tmp3 = Distances(i,j) - dists(i,j);
86 sigma += Weights(i,j)*tmp3*tmp3;
87 }
88 }
89 return sigma;
90 }
91
92 }
93 }
94