Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2015-2023 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 "DimensionalityReductionBase.h"
      23             : #include "reference/ReferenceConfiguration.h"
      24             : #include "core/PlumedMain.h"
      25             : #include "core/Atoms.h"
      26             : 
      27             : namespace PLMD {
      28             : namespace dimred {
      29             : 
      30          44 : void DimensionalityReductionBase::registerKeywords( Keywords& keys ) {
      31          44 :   analysis::AnalysisBase::registerKeywords( keys );
      32          88 :   keys.add("compulsory","NLOW_DIM","number of low-dimensional coordinates required");
      33          88 :   keys.addOutputComponent("coord","default","the low-dimensional projections of the various input configurations");
      34          44 : }
      35             : 
      36          16 : DimensionalityReductionBase::DimensionalityReductionBase( const ActionOptions& ao ):
      37             :   Action(ao),
      38             :   analysis::AnalysisBase(ao),
      39          16 :   dimredbase(NULL) {
      40             :   // Check that some dissimilarity information is available
      41          16 :   if( my_input_data ) {
      42          27 :     if( getName()!="PCA" && !dissimilaritiesWereSet() ) {
      43           0 :       error("dissimilarities have not been calcualted in input actions");
      44             :     }
      45             :     // Now we check if the input was a dimensionality reduction object
      46          15 :     dimredbase = dynamic_cast<DimensionalityReductionBase*>( my_input_data );
      47             :   }
      48             : 
      49             :   // Retrieve the dimension in the low dimensionality space
      50          16 :   nlow=0;
      51          16 :   if( dimredbase ) {
      52           5 :     nlow=dimredbase->nlow;
      53           5 :     log.printf("  projecting in %u dimensional space \n",nlow);
      54          22 :   } else if( keywords.exists("NLOW_DIM") ) {
      55          10 :     parse("NLOW_DIM",nlow);
      56          10 :     if( nlow<1 ) {
      57           0 :       error("dimensionality of low dimensional space must be at least one");
      58             :     }
      59          10 :     log.printf("  projecting in %u dimensional space \n",nlow);
      60             :   }
      61             :   // Now add fake components to the underlying ActionWithValue for the arguments
      62             :   std::string num;
      63          45 :   for(unsigned i=0; i<nlow; ++i) {
      64          29 :     Tools::convert(i+1,num);
      65          29 :     addComponent( "coord-" + num );
      66          58 :     componentIsNotPeriodic( "coord-" + num );
      67             :   }
      68          16 : }
      69             : 
      70           2 : std::vector<Value*> DimensionalityReductionBase::getArgumentList() {
      71             :   std::vector<Value*> arglist( analysis::AnalysisBase::getArgumentList() );
      72           6 :   for(unsigned i=0; i<nlow; ++i) {
      73           4 :     arglist.push_back( getPntrToComponent(i) );
      74             :   }
      75           2 :   return arglist;
      76             : }
      77             : 
      78        1050 : void DimensionalityReductionBase::getProjection( const unsigned& idata, std::vector<double>& point, double& weight ) {
      79        1050 :   if( point.size()!=nlow ) {
      80           0 :     point.resize( nlow );
      81             :   }
      82        1050 :   weight = getWeight(idata);
      83        3150 :   for(unsigned i=0; i<nlow; ++i) {
      84        2100 :     point[i]=projections(idata,i);
      85             :   }
      86        1050 : }
      87             : 
      88          19 : void DimensionalityReductionBase::performAnalysis() {
      89          19 :   log.printf("Generating projections required by action %s \n",getLabel().c_str() );
      90             :   // Resize the tempory array (this is used for out of sample)
      91          19 :   dtargets.resize( getNumberOfDataPoints() );
      92             :   // Resize the projections array
      93          19 :   projections.resize( getNumberOfDataPoints(), nlow );
      94             :   // Retrieve the projections from the previous calculation
      95          19 :   if( dimredbase ) {
      96           6 :     std::vector<double> newp( nlow );
      97             :     double w;
      98        1056 :     for(unsigned i=0; i<getNumberOfDataPoints(); ++i) {
      99        1050 :       dimredbase->getProjection( i, newp, w );
     100             :       plumed_dbg_assert( newp.size()==nlow );
     101        3150 :       for(unsigned j=0; j<nlow; ++j) {
     102        2100 :         projections(i,j)=newp[j];
     103             :       }
     104             :     }
     105             :   }
     106             :   // Calculate matrix of dissimilarities
     107          19 :   Matrix<double> targets( getNumberOfDataPoints(), getNumberOfDataPoints() );
     108             :   targets=0;
     109        2686 :   for(unsigned i=1; i<getNumberOfDataPoints(); ++i) {
     110      367354 :     for(unsigned j=0; j<i; ++j) {
     111      364687 :       targets(i,j)=targets(j,i)=getDissimilarity( i, j );
     112             :     }
     113             :   }
     114             :   // This calculates the projections of the points
     115          19 :   calculateProjections( targets, projections );
     116             :   // Now set the projection values in the underlying object
     117          19 :   if( my_input_data ) {
     118        2620 :     for(unsigned idat=0; idat<getNumberOfDataPoints(); ++idat) {
     119        2602 :       analysis::DataCollectionObject& myref=AnalysisBase::getStoredData(idat,false);
     120             :       std::string num;
     121        7804 :       for(unsigned i=0; i<nlow; ++i) {
     122        5202 :         Tools::convert(i+1,num);
     123       10404 :         myref.setArgument( getLabel() + ".coord-" + num, projections(idat,i) );
     124             :       }
     125             :     }
     126             :   }
     127          19 :   log.printf("Generated projections required by action %s \n",getLabel().c_str() );
     128          19 : }
     129             : 
     130           0 : double DimensionalityReductionBase::calculateStress( const std::vector<double>& p, std::vector<double>& d ) {
     131             : 
     132             :   // Zero derivative and stress accumulators
     133           0 :   for(unsigned i=0; i<p.size(); ++i) {
     134           0 :     d[i]=0.0;
     135             :   }
     136             :   double stress=0;
     137           0 :   std::vector<double> dtmp( p.size() );
     138             : 
     139             :   // Now accumulate total stress on system
     140           0 :   for(unsigned i=0; i<dtargets.size(); ++i) {
     141           0 :     if( dtargets[i]<epsilon ) {
     142           0 :       continue ;
     143             :     }
     144             : 
     145             :     // Calculate distance in low dimensional space
     146             :     double dd=0;
     147           0 :     for(unsigned j=0; j<p.size(); ++j) {
     148           0 :       dtmp[j]=p[j]-projections(i,j);
     149           0 :       dd+=dtmp[j]*dtmp[j];
     150             :     }
     151           0 :     dd = std::sqrt(dd);
     152             : 
     153             :     // Now do transformations and calculate differences
     154           0 :     double ddiff = dd - dtargets[i];
     155             : 
     156             :     // Calculate derivatives
     157           0 :     double pref = 2.*getWeight(i) / dd;
     158           0 :     for(unsigned j=0; j<p.size(); ++j) {
     159           0 :       d[j] += pref*ddiff*dtmp[j];
     160             :     }
     161             : 
     162             :     // Accumulate the total stress
     163           0 :     stress += getWeight(i)*ddiff*ddiff;
     164             :   }
     165           0 :   return stress;
     166             : }
     167             : 
     168             : }
     169             : }