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 "SMACOF.h"
      23             : 
      24             : namespace PLMD {
      25             : namespace dimred {
      26             : 
      27           1 : SMACOF::SMACOF( const Value* target)
      28             : {
      29           1 :   std::vector<unsigned> shape( target->getShape() );
      30           1 :   Distances.resize( shape[0], shape[1] ); Weights.resize( shape[0], shape[1] );
      31         501 :   for(unsigned i=0; i<shape[0]; ++i) {
      32      250500 :     for(unsigned j=0; j<shape[1]; ++j) Distances(i,j) = sqrt( target->get(shape[0]*i+j) );
      33             :   }
      34           1 : }
      35             : 
      36           1 : void SMACOF::optimize( const double& tol, const unsigned& maxloops, std::vector<double>& proj ) {
      37           1 :   unsigned M = Distances.nrows(); unsigned nlow=proj.size() / M; Matrix<double> Z( M, nlow );
      38             :   // Transfer initial projection to matrix
      39             :   unsigned k = 0;
      40         501 :   for(unsigned i=0; i<M; ++i) {
      41        1500 :     for(unsigned j=0; j<nlow; ++j) { Z(i,j)=proj[k]; k++; }
      42             :   }
      43             : 
      44             :   // Calculate V
      45             :   Matrix<double> V(M,M);
      46         501 :   for(unsigned i=0; i<M; ++i) {
      47      250500 :     for(unsigned j=0; j<M; ++j) {
      48      250000 :       if(i==j) continue;
      49      249500 :       V(i,j)=-Weights(i,j);
      50             :     }
      51      250500 :     for(unsigned j=0; j<M; ++j) {
      52      250000 :       if(i==j)continue;
      53      249500 :       V(i,i)-=V(i,j);
      54             :     }
      55             :   }
      56             :   // And pseudo invert V
      57           1 :   Matrix<double> mypseudo(M, M); pseudoInvert(V, mypseudo);
      58           1 :   Matrix<double> dists( M, M ); double myfirstsig = calculateSigma( Z, dists );
      59             : 
      60             :   // initial sigma is made up of the original distances minus the distances between the projections all squared.
      61             :   Matrix<double> temp( M, M ), BZ( M, M ), newZ( M, nlow );
      62          14 :   for(unsigned n=0; n<maxloops; ++n) {
      63          14 :     if(n==maxloops-1) plumed_merror("ran out of steps in SMACOF algorithm");
      64             : 
      65             :     // Recompute BZ matrix
      66        7014 :     for(unsigned i=0; i<M; ++i) {
      67     3507000 :       for(unsigned j=0; j<M; ++j) {
      68     3500000 :         if(i==j) continue;  //skips over the diagonal elements
      69             : 
      70     3493000 :         if( dists(i,j)>0 ) BZ(i,j) = -Weights(i,j)*Distances(i,j) / dists(i,j);
      71           0 :         else BZ(i,j)=0.;
      72             :       }
      73             :       //the diagonal elements are -off diagonal elements BZ(i,i)-=BZ(i,j)   (Equation 8.25)
      74        7000 :       BZ(i,i)=0; //create the space memory for the diagonal elements which are scalars
      75     3507000 :       for(unsigned j=0; j<M; ++j) {
      76     3500000 :         if(i==j) continue;
      77     3493000 :         BZ(i,i)-=BZ(i,j);
      78             :       }
      79             :     }
      80             : 
      81          14 :     mult( mypseudo, BZ, temp); mult(temp, Z, newZ);
      82             :     //Compute new sigma
      83          14 :     double newsig = calculateSigma( newZ, dists );
      84             :     //Computing whether the algorithm has converged (has the mass of the potato changed
      85             :     //when we put it back in the oven!)
      86          14 :     if( fabs( newsig - myfirstsig )<tol ) break;
      87             :     myfirstsig=newsig; Z = newZ;
      88             :   }
      89             : 
      90             :   // Transfer final projection matrix to output proj
      91             :   k = 0;
      92         501 :   for(unsigned i=0; i<M; ++i) {
      93        1500 :     for(unsigned j=0; j<nlow; ++j) { proj[k]=Z(i,j); k++; }
      94             :   }
      95           1 : }
      96             : 
      97          15 : double SMACOF::calculateSigma( const Matrix<double>& Z, Matrix<double>& dists ) {
      98             :   unsigned M = Distances.nrows(); double sigma=0; double totalWeight=0;
      99        7500 :   for(unsigned i=1; i<M; ++i) {
     100     1878735 :     for(unsigned j=0; j<i; ++j) {
     101     5613750 :       double dlow=0; for(unsigned k=0; k<Z.ncols(); ++k) { double tmp=Z(i,k) - Z(j,k); dlow+=tmp*tmp; }
     102     1871250 :       dists(i,j)=dists(j,i)=sqrt(dlow); double tmp3 = Distances(i,j) - dists(i,j);
     103     1871250 :       sigma += Weights(i,j)*tmp3*tmp3; totalWeight+=Weights(i,j);
     104             :     }
     105             :   }
     106          15 :   return sigma / totalWeight;
     107             : }
     108             : 
     109             : }
     110             : }