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 "matrixtools/MatrixOperationBase.h"
      23             : #include "core/ActionRegister.h"
      24             : #include "tools/Random.h"
      25             : 
      26             : //+PLUMEDOC LANDMARKS FARTHEST_POINT_SAMPLING
      27             : /*
      28             : Select a set of landmarks using farthest point sampling.
      29             : 
      30             : \par Examples
      31             : 
      32             : */
      33             : //+ENDPLUMEDOC
      34             : 
      35             : namespace PLMD {
      36             : namespace landmarks {
      37             : 
      38             : class FarthestPointSampling : public matrixtools::MatrixOperationBase {
      39             : private:
      40             :   unsigned seed;
      41             :   unsigned nlandmarks;
      42             : public:
      43             :   static void registerKeywords( Keywords& keys );
      44             :   explicit FarthestPointSampling( const ActionOptions& ao );
      45           0 :   unsigned getNumberOfDerivatives() override {
      46           0 :     return 0;
      47             :   }
      48             :   void prepare() override ;
      49             :   void calculate() override ;
      50           0 :   void apply() override {}
      51           0 :   double getForceOnMatrixElement( const unsigned& jrow, const unsigned& krow ) const {
      52           0 :     plumed_merror("this should not be called");
      53             :   }
      54             : };
      55             : 
      56             : PLUMED_REGISTER_ACTION(FarthestPointSampling,"FARTHEST_POINT_SAMPLING")
      57             : 
      58           4 : void FarthestPointSampling::registerKeywords( Keywords& keys ) {
      59           4 :   matrixtools::MatrixOperationBase::registerKeywords( keys );
      60           8 :   keys.add("compulsory","NZEROS","the number of landmark points that you want to select");
      61           8 :   keys.add("compulsory","SEED","1234","a random number seed");
      62           4 :   keys.setValueDescription("a vector which has as many elements as there are rows in the input matrix of dissimilarities. NZEROS of the elements in this vector are equal to one, the rest of the elements are equal to zero.  The nodes that have elements equal to one are the NZEROS points that are farthest appart according to the input dissimilarities");
      63           4 : }
      64             : 
      65           1 : FarthestPointSampling::FarthestPointSampling( const ActionOptions& ao ):
      66             :   Action(ao),
      67           1 :   MatrixOperationBase(ao) {
      68           1 :   if( getPntrToArgument(0)->getShape()[0]!=getPntrToArgument(0)->getShape()[1] ) {
      69           0 :     error("input to this argument should be a square matrix of dissimilarities");
      70             :   }
      71           1 :   parse("NZEROS",nlandmarks);
      72           1 :   parse("SEED",seed);
      73           1 :   log.printf("  selecting %d landmark points \n", nlandmarks );
      74             : 
      75           1 :   std::vector<unsigned> shape(1);
      76           1 :   shape[0] = getPntrToArgument(0)->getShape()[0];
      77           1 :   addValue( shape );
      78           1 :   setNotPeriodic();
      79           1 :   getPntrToComponent(0)->buildDataStore();
      80           1 : }
      81             : 
      82           1 : void FarthestPointSampling::prepare() {
      83           1 :   Value* myval = getPntrToComponent(0);
      84           1 :   if( myval->getShape()[0]!=getPntrToArgument(0)->getShape()[0] ) {
      85           1 :     std::vector<unsigned> shape(1);
      86           1 :     shape[0] = getPntrToArgument(0)->getShape()[0];
      87           1 :     myval->setShape(shape);
      88             :   }
      89           4 :   for(unsigned i=0; i<nlandmarks; ++i) {
      90           3 :     myval->set( i, 0.0 );
      91             :   }
      92          11 :   for(unsigned i=nlandmarks; i<myval->getShape()[0]; ++i) {
      93          10 :     myval->set( i, 1.0 );
      94             :   }
      95           1 : }
      96             : 
      97           1 : void FarthestPointSampling::calculate() {
      98           1 :   Value* myval=getPntrToComponent(0);
      99           1 :   unsigned npoints = getPntrToArgument(0)->getShape()[0];
     100          14 :   for(unsigned i=0; i<npoints; ++i) {
     101          13 :     myval->set( i, 1.0 );
     102             :   }
     103           1 :   std::vector<unsigned> landmarks( nlandmarks );
     104             : 
     105             :   // Select first point at random
     106           1 :   Random random;
     107           1 :   random.setSeed(-seed);
     108           1 :   double rand=random.RandU01();
     109           1 :   landmarks[0] = std::floor( npoints*rand );
     110           1 :   myval->set( landmarks[0], 0 );
     111             : 
     112             :   // Now find distance to all other points (N.B. We can use squared distances here for speed)
     113           1 :   Matrix<double> distances( nlandmarks, npoints );
     114             :   Value* myarg = getPntrToArgument(0);
     115          14 :   for(unsigned i=0; i<npoints; ++i) {
     116          13 :     distances(0,i) = myarg->get( landmarks[0]*npoints + i );
     117             :   }
     118             : 
     119             :   // Now find all other landmarks
     120           3 :   for(unsigned i=1; i<nlandmarks; ++i) {
     121             :     // Find point that has the largest minimum distance from the landmarks selected thus far
     122             :     double maxd=0;
     123          28 :     for(unsigned j=0; j<npoints; ++j) {
     124          26 :       double mind=distances(0,j);
     125          39 :       for(unsigned k=1; k<i; ++k) {
     126          13 :         if( distances(k,j)<mind ) {
     127             :           mind=distances(k,j);
     128             :         }
     129             :       }
     130          26 :       if( mind>maxd ) {
     131             :         maxd=mind;
     132           5 :         landmarks[i]=j;
     133             :       }
     134             :     }
     135           2 :     myval->set( landmarks[i], 0 );
     136          28 :     for(unsigned k=0; k<npoints; ++k) {
     137          26 :       distances(i,k) = myarg->get( landmarks[i]*npoints + k );
     138             :     }
     139             :   }
     140           1 : }
     141             : 
     142             : }
     143             : }