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 "core/ActionRegister.h"
      23             : #include "core/PlumedMain.h"
      24             : #include "core/ActionSet.h"
      25             : #include "tools/Random.h"
      26             : #include "tools/ConjugateGradient.h"
      27             : #include "analysis/AnalysisBase.h"
      28             : #include "reference/ReferenceConfiguration.h"
      29             : #include "DimensionalityReductionBase.h"
      30             : #include "PCA.h"
      31             : 
      32             : //+PLUMEDOC DIMRED PROJECT_ALL_ANALYSIS_DATA
      33             : /*
      34             : Find projections of all non-landmark points using the embedding calculated by a dimensionality reduction optimization calculation.
      35             : 
      36             : \par Examples
      37             : 
      38             : */
      39             : //+ENDPLUMEDOC
      40             : 
      41             : namespace PLMD {
      42             : namespace dimred {
      43             : 
      44             : class ProjectNonLandmarkPoints : public analysis::AnalysisBase {
      45             : private:
      46             : /// Tolerance for conjugate gradient algorithm
      47             :   double cgtol;
      48             : /// Number of diemsions in low dimensional space
      49             :   unsigned nlow;
      50             : /// The class that calculates the projection of the data that is required
      51             :   DimensionalityReductionBase* mybase;
      52             : /// Generate a projection of the ith data point - this is called in two routine
      53             :   void generateProjection( const unsigned& idat, std::vector<double>& point );
      54             : public:
      55             :   static void registerKeywords( Keywords& keys );
      56             :   explicit ProjectNonLandmarkPoints( const ActionOptions& ao );
      57             : /// Get a reference configuration (this returns the projection)
      58             :   analysis::DataCollectionObject& getStoredData( const unsigned& idat, const bool& calcdist );
      59             : /// Overwrite getArguments so we get arguments from underlying class
      60             :   std::vector<Value*> getArgumentList();
      61             : /// This does nothing -- projections are calculated when getDataPoint and getReferenceConfiguration are called
      62           2 :   void performAnalysis() {}
      63             : /// This just calls calculate stress in the underlying projection object
      64             :   double calculateStress( const std::vector<double>& pp, std::vector<double>& der );
      65             : /// Overwrite virtual function in ActionWithVessel
      66           0 :   void performTask( const unsigned&, const unsigned&, MultiValue& ) const {
      67           0 :     plumed_error();
      68             :   }
      69             : };
      70             : 
      71       13789 : PLUMED_REGISTER_ACTION(ProjectNonLandmarkPoints,"PROJECT_ALL_ANALYSIS_DATA")
      72             : 
      73           6 : void ProjectNonLandmarkPoints::registerKeywords( Keywords& keys ) {
      74           6 :   analysis::AnalysisBase::registerKeywords( keys );
      75          12 :   keys.add("compulsory","PROJECTION","the projection that you wish to generate out-of-sample projections with");
      76          12 :   keys.add("compulsory","CGTOL","1E-6","the tolerance for the conjugate gradient optimization");
      77          12 :   keys.addOutputComponent("coord","default","the low-dimensional projections of the various input configurations");
      78           6 : }
      79             : 
      80           2 : ProjectNonLandmarkPoints::ProjectNonLandmarkPoints( const ActionOptions& ao ):
      81             :   Action(ao),
      82             :   analysis::AnalysisBase(ao),
      83           2 :   mybase(NULL) {
      84             :   std::string myproj;
      85           2 :   parse("PROJECTION",myproj);
      86           2 :   mybase = plumed.getActionSet().selectWithLabel<DimensionalityReductionBase*>( myproj );
      87           2 :   if( !mybase ) {
      88           0 :     error("could not find projection of data named " + myproj );
      89             :   }
      90             :   // Add the dependency and set the dimensionality
      91           2 :   addDependency( mybase );
      92           2 :   nlow = mybase->nlow;
      93             :   // Add fake components to the underlying ActionWithValue for the arguments
      94             :   std::string num;
      95           6 :   for(unsigned i=0; i<nlow; ++i) {
      96           4 :     Tools::convert(i+1,num);
      97           4 :     addComponent( "coord-" + num );
      98           8 :     componentIsNotPeriodic( "coord-" + num );
      99             :   }
     100             : 
     101           2 :   log.printf("  generating out-of-sample projections using projection with label %s \n",myproj.c_str() );
     102           4 :   parse("CGTOL",cgtol);
     103           2 : }
     104             : 
     105           0 : std::vector<Value*> ProjectNonLandmarkPoints::getArgumentList() {
     106             :   std::vector<Value*> arglist( analysis::AnalysisBase::getArgumentList() );
     107           0 :   for(unsigned i=0; i<nlow; ++i) {
     108           0 :     arglist.push_back( getPntrToComponent(i) );
     109             :   }
     110           0 :   return arglist;
     111             : }
     112             : 
     113        1046 : void ProjectNonLandmarkPoints::generateProjection( const unsigned& idat, std::vector<double>& point ) {
     114        1046 :   PCA* ispca = dynamic_cast<PCA*>( mybase );
     115        1046 :   if( ispca ) {
     116         546 :     ispca->getProjection( my_input_data->getStoredData(idat,false), point );
     117             :   } else {
     118             :     ConjugateGradient<ProjectNonLandmarkPoints> myminimiser( this );
     119             :     unsigned closest=0;
     120         500 :     double mindist = std::sqrt( getDissimilarity( mybase->getDataPointIndexInBase(0), idat ) );
     121         500 :     mybase->setTargetDistance( 0, mindist );
     122      125000 :     for(unsigned i=1; i<mybase->getNumberOfDataPoints(); ++i) {
     123      124500 :       double dist = std::sqrt( getDissimilarity( mybase->getDataPointIndexInBase(i), idat ) );
     124      124500 :       mybase->setTargetDistance( i, dist );
     125      124500 :       if( dist<mindist ) {
     126        2496 :         mindist=dist;
     127        2496 :         closest=i;
     128             :       }
     129             :     }
     130             :     // Put the initial guess near to the closest landmark  -- may wish to use grid here again Sandip??
     131         500 :     Random random;
     132         500 :     random.setSeed(-1234);
     133        1500 :     for(unsigned j=0; j<nlow; ++j) {
     134        1000 :       point[j]=mybase->projections(closest,j) + (random.RandU01() - 0.5)*0.01;
     135             :     }
     136         500 :     myminimiser.minimise( cgtol, point, &ProjectNonLandmarkPoints::calculateStress );
     137             :   }
     138        1046 : }
     139             : 
     140        1046 : analysis::DataCollectionObject& ProjectNonLandmarkPoints::getStoredData( const unsigned& idat, const bool& calcdist ) {
     141        1046 :   std::vector<double> pp(nlow);
     142        1046 :   generateProjection( idat, pp );
     143             :   std::string num;
     144             :   analysis::DataCollectionObject& myref=AnalysisBase::getStoredData(idat,calcdist);
     145        3138 :   for(unsigned i=0; i<nlow; ++i) {
     146        2092 :     Tools::convert(i+1,num);
     147        4184 :     myref.setArgument( getLabel() + ".coord-" + num, pp[i] );
     148             :   }
     149        1046 :   return myref;
     150             : }
     151             : 
     152       16744 : double ProjectNonLandmarkPoints::calculateStress( const std::vector<double>& pp, std::vector<double>& der ) {
     153       16744 :   return mybase->calculateStress( pp, der );
     154             : }
     155             : 
     156             : }
     157             : }
     158             :