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 { plumed_error(); }
      67             : };
      68             : 
      69       10423 : PLUMED_REGISTER_ACTION(ProjectNonLandmarkPoints,"PROJECT_ALL_ANALYSIS_DATA")
      70             : 
      71           3 : void ProjectNonLandmarkPoints::registerKeywords( Keywords& keys ) {
      72           3 :   analysis::AnalysisBase::registerKeywords( keys );
      73           6 :   keys.add("compulsory","PROJECTION","the projection that you wish to generate out-of-sample projections with");
      74           6 :   keys.add("compulsory","CGTOL","1E-6","the tolerance for the conjugate gradient optimization");
      75           6 :   keys.addOutputComponent("coord","default","the low-dimensional projections of the various input configurations");
      76           3 : }
      77             : 
      78           2 : ProjectNonLandmarkPoints::ProjectNonLandmarkPoints( const ActionOptions& ao ):
      79             :   Action(ao),
      80             :   analysis::AnalysisBase(ao),
      81           2 :   mybase(NULL)
      82             : {
      83           2 :   std::string myproj; parse("PROJECTION",myproj);
      84           2 :   mybase = plumed.getActionSet().selectWithLabel<DimensionalityReductionBase*>( myproj );
      85           2 :   if( !mybase ) error("could not find projection of data named " + myproj );
      86             :   // Add the dependency and set the dimensionality
      87           2 :   addDependency( mybase ); nlow = mybase->nlow;
      88             :   // Add fake components to the underlying ActionWithValue for the arguments
      89             :   std::string num;
      90           6 :   for(unsigned i=0; i<nlow; ++i) {
      91          12 :     Tools::convert(i+1,num); addComponent( "coord-" + num ); componentIsNotPeriodic( "coord-" + num );
      92             :   }
      93             : 
      94           2 :   log.printf("  generating out-of-sample projections using projection with label %s \n",myproj.c_str() );
      95           4 :   parse("CGTOL",cgtol);
      96           2 : }
      97             : 
      98           0 : std::vector<Value*> ProjectNonLandmarkPoints::getArgumentList() {
      99             :   std::vector<Value*> arglist( analysis::AnalysisBase::getArgumentList() );
     100           0 :   for(unsigned i=0; i<nlow; ++i) arglist.push_back( getPntrToComponent(i) );
     101           0 :   return arglist;
     102             : }
     103             : 
     104        1046 : void ProjectNonLandmarkPoints::generateProjection( const unsigned& idat, std::vector<double>& point ) {
     105        1046 :   PCA* ispca = dynamic_cast<PCA*>( mybase );
     106        1046 :   if( ispca ) {
     107         546 :     ispca->getProjection( my_input_data->getStoredData(idat,false), point );
     108             :   } else {
     109             :     ConjugateGradient<ProjectNonLandmarkPoints> myminimiser( this );
     110         500 :     unsigned closest=0; double mindist = std::sqrt( getDissimilarity( mybase->getDataPointIndexInBase(0), idat ) );
     111         500 :     mybase->setTargetDistance( 0, mindist );
     112      125000 :     for(unsigned i=1; i<mybase->getNumberOfDataPoints(); ++i) {
     113      124500 :       double dist = std::sqrt( getDissimilarity( mybase->getDataPointIndexInBase(i), idat ) );
     114      124500 :       mybase->setTargetDistance( i, dist );
     115      124500 :       if( dist<mindist ) { mindist=dist; closest=i; }
     116             :     }
     117             :     // Put the initial guess near to the closest landmark  -- may wish to use grid here again Sandip??
     118         500 :     Random random; random.setSeed(-1234);
     119        1500 :     for(unsigned j=0; j<nlow; ++j) point[j]=mybase->projections(closest,j) + (random.RandU01() - 0.5)*0.01;
     120         500 :     myminimiser.minimise( cgtol, point, &ProjectNonLandmarkPoints::calculateStress );
     121             :   }
     122        1046 : }
     123             : 
     124        1046 : analysis::DataCollectionObject& ProjectNonLandmarkPoints::getStoredData( const unsigned& idat, const bool& calcdist ) {
     125        1046 :   std::vector<double> pp(nlow); generateProjection( idat, pp ); std::string num;
     126             :   analysis::DataCollectionObject& myref=AnalysisBase::getStoredData(idat,calcdist);
     127        3138 :   for(unsigned i=0; i<nlow; ++i) { Tools::convert(i+1,num); myref.setArgument( getLabel() + ".coord-" + num, pp[i] ); }
     128        1046 :   return myref;
     129             : }
     130             : 
     131       16744 : double ProjectNonLandmarkPoints::calculateStress( const std::vector<double>& pp, std::vector<double>& der ) {
     132       16744 :   return mybase->calculateStress( pp, der );
     133             : }
     134             : 
     135             : }
     136             : }
     137             :