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