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 { return 0; }
      46             :   void prepare() override ;
      47             :   void calculate() override ;
      48           0 :   void apply() override {}
      49           0 :   double getForceOnMatrixElement( const unsigned& jrow, const unsigned& krow ) const { plumed_merror("this should not be called"); }
      50             : };
      51             : 
      52             : PLUMED_REGISTER_ACTION(FarthestPointSampling,"FARTHEST_POINT_SAMPLING")
      53             : 
      54           4 : void FarthestPointSampling::registerKeywords( Keywords& keys ) {
      55           4 :   matrixtools::MatrixOperationBase::registerKeywords( keys );
      56           8 :   keys.add("compulsory","NZEROS","the number of landmark points that you want to select");
      57           8 :   keys.add("compulsory","SEED","1234","a random number seed");
      58           8 :   keys.setValueDescription("vector","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");
      59           4 : }
      60             : 
      61           1 : FarthestPointSampling::FarthestPointSampling( const ActionOptions& ao ):
      62             :   Action(ao),
      63           1 :   MatrixOperationBase(ao)
      64             : {
      65           1 :   if( getPntrToArgument(0)->getShape()[0]!=getPntrToArgument(0)->getShape()[1] ) error("input to this argument should be a square matrix of dissimilarities");
      66           2 :   parse("NZEROS",nlandmarks); parse("SEED",seed);
      67           1 :   log.printf("  selecting %d landmark points \n", nlandmarks );
      68             : 
      69           1 :   std::vector<unsigned> shape(1); shape[0] = getPntrToArgument(0)->getShape()[0];
      70           1 :   addValue( shape ); setNotPeriodic(); getPntrToComponent(0)->buildDataStore();
      71           1 : }
      72             : 
      73           1 : void FarthestPointSampling::prepare() {
      74           1 :   Value* myval = getPntrToComponent(0);
      75           1 :   if( myval->getShape()[0]!=getPntrToArgument(0)->getShape()[0] ) {
      76           1 :     std::vector<unsigned> shape(1); shape[0] = getPntrToArgument(0)->getShape()[0]; myval->setShape(shape);
      77             :   }
      78           4 :   for(unsigned i=0; i<nlandmarks; ++i) myval->set( i, 0.0 );
      79          11 :   for(unsigned i=nlandmarks; i<myval->getShape()[0]; ++i) myval->set( i, 1.0 );
      80           1 : }
      81             : 
      82           1 : void FarthestPointSampling::calculate() {
      83           1 :   Value* myval=getPntrToComponent(0); unsigned npoints = getPntrToArgument(0)->getShape()[0];
      84          14 :   for(unsigned i=0; i<npoints; ++i) myval->set( i, 1.0 );
      85           1 :   std::vector<unsigned> landmarks( nlandmarks );
      86             : 
      87             :   // Select first point at random
      88           1 :   Random random; random.setSeed(-seed); double rand=random.RandU01();
      89           1 :   landmarks[0] = std::floor( npoints*rand ); myval->set( landmarks[0], 0 );
      90             : 
      91             :   // Now find distance to all other points (N.B. We can use squared distances here for speed)
      92           1 :   Matrix<double> distances( nlandmarks, npoints ); Value* myarg = getPntrToArgument(0);
      93          14 :   for(unsigned i=0; i<npoints; ++i) distances(0,i) = myarg->get( landmarks[0]*npoints + i );
      94             : 
      95             :   // Now find all other landmarks
      96           3 :   for(unsigned i=1; i<nlandmarks; ++i) {
      97             :     // Find point that has the largest minimum distance from the landmarks selected thus far
      98             :     double maxd=0;
      99          28 :     for(unsigned j=0; j<npoints; ++j) {
     100          26 :       double mind=distances(0,j);
     101          39 :       for(unsigned k=1; k<i; ++k) {
     102          13 :         if( distances(k,j)<mind ) { mind=distances(k,j); }
     103             :       }
     104          26 :       if( mind>maxd ) { maxd=mind; landmarks[i]=j; }
     105             :     }
     106           2 :     myval->set( landmarks[i], 0 );
     107          28 :     for(unsigned k=0; k<npoints; ++k) distances(i,k) = myarg->get( landmarks[i]*npoints + k );
     108             :   }
     109           1 : }
     110             : 
     111             : }
     112             : }