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