Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2016-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 "core/ActionWithVector.h"
      23             : #include "core/ActionRegister.h"
      24             : #include "tools/ConjugateGradient.h"
      25             : #include "tools/SwitchingFunction.h"
      26             : #include "tools/OpenMP.h"
      27             : #include "tools/Random.h"
      28             : 
      29             : namespace PLMD {
      30             : namespace dimred {
      31             : 
      32             : //+PLUMEDOC DIMRED PROJECT_POINTS
      33             : /*
      34             : Find the projection of a point in a low dimensional space by matching the (transformed) distance between it and a series of reference configurations that were input
      35             : 
      36             : \par Examples
      37             : 
      38             : */
      39             : //+ENDPLUMEDOC
      40             : 
      41             : class ProjectPoints : public ActionWithVector {
      42             : private:
      43             :   double cgtol;
      44             :   unsigned dimout;
      45             :   mutable std::vector<unsigned> rowstart;
      46             :   std::vector<SwitchingFunction> switchingFunction;
      47             :   ConjugateGradient<ProjectPoints> myminimiser;
      48             :   void getProjection( const unsigned& current, std::vector<double>& point ) const ;
      49             : public:
      50             :   static void registerKeywords( Keywords& keys );
      51             :   ProjectPoints( const ActionOptions& );
      52           0 :   unsigned getNumberOfDerivatives() { return 0; }
      53             :   void prepare() override ;
      54             :   void performTask( const unsigned& current, MultiValue& myvals ) const override ;
      55             :   double calculateStress( const std::vector<double>& pp, std::vector<double>& der );
      56             :   void calculate() override ;
      57         168 :   void apply() override {}
      58             : };
      59             : 
      60             : PLUMED_REGISTER_ACTION(ProjectPoints,"PROJECT_POINTS")
      61             : 
      62           6 : void ProjectPoints::registerKeywords( Keywords& keys ) {
      63           6 :   ActionWithVector::registerKeywords( keys ); keys.use("ARG");
      64          12 :   keys.add("numbered","TARGET","the matrix of target quantities that you would like to match");
      65          12 :   keys.add("numbered","FUNC","a function that is applied on the distances between the points in the low dimensional space");
      66          12 :   keys.add("numbered","WEIGHTS","the matrix with the weights of the target quantities");
      67          12 :   keys.add("compulsory","CGTOL","1E-6","the tolerance for the conjugate gradient minimization");
      68          12 :   keys.addOutputComponent("coord","default","the coordinates of the points in the low dimensional space");
      69           6 : }
      70             : 
      71             : 
      72           2 : ProjectPoints::ProjectPoints( const ActionOptions& ao ) :
      73             :   Action(ao),
      74             :   ActionWithVector(ao),
      75           2 :   rowstart(OpenMP::getNumThreads()),
      76           4 :   myminimiser( this )
      77             : {
      78           2 :   dimout = getNumberOfArguments(); unsigned nvals=getPntrToArgument(0)->getNumberOfValues();
      79           6 :   for(unsigned i=0; i<getNumberOfArguments(); ++i) {
      80           4 :     if( nvals!=getPntrToArgument(i)->getNumberOfValues() ) error("mismatch between numbers of projections");
      81             :   }
      82           2 :   std::vector<Value*> args( getArguments() ), target, weights; std::string sfd, errors; unsigned ntoproj=0;
      83             :   // Read in target "distances" and target weights
      84           2 :   for(unsigned i=1;; ++i) {
      85           8 :     target.resize(0); if( !parseArgumentList("TARGET",i,target) ) break;
      86           2 :     std::string inum; Tools::convert( i, inum );
      87           2 :     if( target.size()!=1 ) error("should only be one value in input to TARGET" + inum );
      88           2 :     if( (target[0]->getRank()!=1 && target[0]->getRank()!=2) || target[0]->hasDerivatives() ) error("input to TARGET" + inum + " keyword should be a vector/matrix");
      89           2 :     if( target[0]->getShape()[0]!=nvals ) error("number of rows in target matrix should match number of input coordinates");
      90           2 :     if( i==1 && target[0]->getRank()==1 ) { ntoproj = 1; }
      91           1 :     else if( ntoproj==1 && target[0]->getRank()!=1 ) error("mismatch between numbers of target distances");
      92           1 :     else if( i==1 ) ntoproj = target[0]->getShape()[1];
      93           0 :     else if( ntoproj!=target[0]->getShape()[1] ) error("mismatch between numbers of target distances");
      94           4 :     if( !parseArgumentList("WEIGHTS",i,weights) ) error("missing WEIGHTS" + inum + " keyword in input");
      95           2 :     if( weights.size()!=1 ) error("should only be one value in input to WEIGHTS" + inum );
      96           2 :     if( weights[0]->getRank()!=1 || weights[0]->hasDerivatives() ) error("input to WEIGHTS" + inum + " keyword should be a vector");
      97           2 :     if( weights[0]->getShape()[0]!=nvals ) error("number of weights should match number of input coordinates");
      98           2 :     target[0]->buildDataStore(); weights[0]->buildDataStore(); args.push_back( target[0] ); args.push_back( weights[0] );
      99           4 :     bool has_sf = parseNumbered("FUNC",i,sfd); switchingFunction.push_back( SwitchingFunction() );
     100           2 :     if( !has_sf ) {
     101           0 :       switchingFunction[i-1].set( "CUSTOM FUNC=1-sqrt(x2) R_0=1.0", errors );
     102             :     } else {
     103           2 :       switchingFunction[i-1].set( sfd, errors );
     104           2 :       if( errors.length()!=0 ) error("problem reading switching function description " + errors);
     105             :     }
     106           2 :     log.printf("  %sth term seeks to match tranformed distances with those in matrix %s \n", inum.c_str(), target[0]->getName().c_str() );
     107           2 :     log.printf("  in %sth term distances are transformed by 1-switching function with r_0=%s \n", inum.c_str(), switchingFunction[i-1].description().c_str() );
     108           2 :     log.printf("  in %sth term weights of matrix elements in stress function are given by %s \n", inum.c_str(), weights[0]->getName().c_str() );
     109           2 :   }
     110           2 :   std::vector<unsigned> shape(1); shape[0]=ntoproj; if( ntoproj==1 ) shape.resize(0);
     111           6 :   for(unsigned i=0; i<dimout; ++i) {
     112           4 :     std::string num; Tools::convert( i+1, num ); addComponent( "coord-" + num, shape );
     113           8 :     componentIsNotPeriodic( "coord-" + num );
     114             :   }
     115             :   // Create a list of tasks to perform
     116           4 :   parse("CGTOL",cgtol); log.printf("  tolerance for conjugate gradient algorithm equals %f \n",cgtol);
     117           2 :   requestArguments( args ); checkRead();
     118           2 : }
     119             : 
     120         169 : void ProjectPoints::prepare() {
     121         169 :   if( getPntrToComponent(0)->getRank()==0 ) return;
     122             : 
     123           1 :   std::vector<unsigned> shape(1); shape[0] = getPntrToArgument(dimout)->getShape()[0];
     124           3 :   for(unsigned i=0; i<dimout; ++i) {
     125           2 :     if( getPntrToComponent(i)->getShape()[0]!=shape[0] ) getPntrToComponent(i)->setShape(shape);
     126             :   }
     127             : }
     128             : 
     129       24180 : double ProjectPoints::calculateStress( const std::vector<double>& pp, std::vector<double>& der ) {
     130       24180 :   unsigned nmatrices = ( getNumberOfArguments() - dimout ) / 2; double stress=0;
     131             : 
     132       24180 :   unsigned t=OpenMP::getThreadNum();
     133       24180 :   std::vector<double> dtmp( pp.size() ); unsigned nland = getPntrToArgument(0)->getShape()[0];
     134     4889418 :   for(unsigned i=0; i<nland; ++i) {
     135             :     // Calculate distance in low dimensional space
     136    14595714 :     double dd2 = 0; for(unsigned k=0; k<pp.size(); ++k) { dtmp[k] = pp[k] - getPntrToArgument(k)->get(i); dd2 += dtmp[k]*dtmp[k]; }
     137             : 
     138     9730476 :     for(unsigned k=0; k<nmatrices; ++k ) {
     139             :       // Now do transformations and calculate differences
     140     4865238 :       double df, fd = 1. - switchingFunction[k].calculateSqr( dd2, df );
     141             :       // Get the weight for this connection
     142     4865238 :       double weight = getPntrToArgument( dimout + 2*k + 1 )->get( i );
     143             :       // Get the difference for the connection
     144     4865238 :       double fdiff = fd - getPntrToArgument( dimout + 2*k )->get( rowstart[t]+i );
     145             :       // Calculate derivatives
     146    14595714 :       double pref = -2.*weight*fdiff*df; for(unsigned n=0; n<pp.size(); ++n) der[n]+=pref*dtmp[n];
     147             :       // Accumulate the total stress
     148     4865238 :       stress += weight*fdiff*fdiff;
     149             :     }
     150             :   }
     151       24180 :   return stress;
     152             : }
     153             : 
     154         668 : void ProjectPoints::getProjection( const unsigned& current, std::vector<double>& point ) const {
     155         668 :   Value* targ = getPntrToArgument( dimout ); unsigned nland = getPntrToArgument(0)->getShape()[0];
     156         668 :   unsigned base = current; if( targ->getRank()==2 ) base = current*targ->getShape()[1];
     157         668 :   unsigned closest=0; double mindist = targ->get( base );
     158      139112 :   for(unsigned i=1; i<nland; ++i) {
     159      138444 :     double dist = targ->get( base + i );
     160      138444 :     if( dist<mindist ) { mindist=dist; closest=i; }
     161             :   }
     162             :   // Put the initial guess near to the closest landmark  -- may wish to use grid here again Sandip??
     163         668 :   Random random; random.setSeed(-1234);
     164        2004 :   for(unsigned j=0; j<dimout; ++j) point[j] = getPntrToArgument(j)->get(closest) + (random.RandU01() - 0.5)*0.01;
     165             :   // And do the optimisation
     166         668 :   rowstart[OpenMP::getThreadNum()]=current; if( targ->getRank()==2 ) rowstart[OpenMP::getThreadNum()] = current*targ->getShape()[1];
     167         668 :   myminimiser.minimise( cgtol, point, &ProjectPoints::calculateStress );
     168         668 : }
     169             : 
     170         500 : void ProjectPoints::performTask( const unsigned& current, MultiValue& myvals ) const {
     171         500 :   std::vector<double> point( dimout ); getProjection( current, point );
     172        1500 :   for(unsigned j=0; j<dimout; ++j) myvals.setValue( getConstPntrToComponent(j)->getPositionInStream(), point[j] );
     173         500 : }
     174             : 
     175         169 : void ProjectPoints::calculate() {
     176         169 :   if( getPntrToComponent(0)->getRank()==0 ) {
     177         168 :     std::vector<double> point( dimout ); getProjection( 0, point );
     178         504 :     for(unsigned i=0; i<dimout; ++i) getPntrToComponent(i)->set(point[i]);
     179           1 :   } else runAllTasks();
     180         169 : }
     181             : 
     182             : }
     183             : }