Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2014-2017 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/ActionWithMatrix.h"
      23             : #include "core/ActionRegister.h"
      24             : 
      25             : //+PLUMEDOC MCOLVAR VSTACK
      26             : /*
      27             : Create a matrix by stacking vectors together
      28             : 
      29             : \par Examples
      30             : 
      31             : */
      32             : //+ENDPLUMEDOC
      33             : 
      34             : namespace PLMD {
      35             : namespace valtools {
      36             : 
      37             : class VStack : public ActionWithMatrix {
      38             : private:
      39             :   std::vector<bool> stored;
      40             : public:
      41             :   static void registerKeywords( Keywords& keys );
      42             : /// Constructor
      43             :   explicit VStack(const ActionOptions&);
      44             : /// Get the number of derivatives
      45         294 :   unsigned getNumberOfDerivatives() override {
      46         294 :     return 0;
      47             :   }
      48             : ///
      49             :   void prepare() override ;
      50             : ///
      51     1270505 :   unsigned getNumberOfColumns() const override {
      52     1270505 :     return getNumberOfArguments();
      53             :   }
      54             : ///
      55             :   void setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const override ;
      56             : ///
      57             :   void performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const override ;
      58             : ///
      59             :   void runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const override ;
      60             : ///
      61             :   void getMatrixColumnTitles( std::vector<std::string>& argnames ) const override ;
      62             : };
      63             : 
      64             : PLUMED_REGISTER_ACTION(VStack,"VSTACK")
      65             : 
      66         254 : void VStack::registerKeywords( Keywords& keys ) {
      67         254 :   ActionWithMatrix::registerKeywords( keys );
      68         254 :   keys.use("ARG");
      69         254 :   keys.setValueDescription("a matrix that contains the input vectors in its columns");
      70         254 : }
      71             : 
      72         139 : VStack::VStack(const ActionOptions& ao):
      73             :   Action(ao),
      74         139 :   ActionWithMatrix(ao) {
      75         139 :   if( getNumberOfArguments()==0 ) {
      76           0 :     error("no arguments were specificed");
      77             :   }
      78         139 :   if( getPntrToArgument(0)->getRank()>1 ) {
      79           0 :     error("all arguments should be vectors");
      80             :   }
      81             :   unsigned nvals=1;
      82             :   bool periodic=false;
      83             :   std::string smin, smax;
      84         139 :   if( getPntrToArgument(0)->getRank()==1 ) {
      85         119 :     nvals = getPntrToArgument(0)->getShape()[0];
      86             :   }
      87         139 :   if( getPntrToArgument(0)->isPeriodic() ) {
      88             :     periodic=true;
      89          24 :     getPntrToArgument(0)->getDomain( smin, smax );
      90             :   }
      91             : 
      92        1228 :   for(unsigned i=0; i<getNumberOfArguments(); ++i) {
      93        1089 :     if( getPntrToArgument(i)->getRank()>1 || (getPntrToArgument(i)->getRank()==1 && getPntrToArgument(i)->hasDerivatives()) ) {
      94           0 :       error("all arguments should be vectors");
      95             :     }
      96        1089 :     if( getPntrToArgument(i)->getRank()==0 ) {
      97          41 :       if( nvals!=1 ) {
      98           0 :         error("all input vector should have same number of elements");
      99             :       }
     100        1048 :     } else if( getPntrToArgument(i)->getShape()[0]!=nvals ) {
     101           0 :       error("all input vector should have same number of elements");
     102             :     }
     103        1089 :     if( periodic ) {
     104          51 :       if( !getPntrToArgument(i)->isPeriodic() ) {
     105           0 :         error("one argument is periodic but " + getPntrToArgument(i)->getName() + " is not periodic");
     106             :       }
     107             :       std::string tmin, tmax;
     108          51 :       getPntrToArgument(i)->getDomain( tmin, tmax );
     109          51 :       if( tmin!=smin || tmax!=smax ) {
     110           0 :         error("domain of argument " + getPntrToArgument(i)->getName() + " is different from domain for all other arguments");
     111             :       }
     112        1038 :     } else if( getPntrToArgument(i)->isPeriodic() ) {
     113           0 :       error("one argument is not periodic but " + getPntrToArgument(i)->getName() + " is periodic");
     114             :     }
     115             :   }
     116             :   // And create a value to hold the matrix
     117         139 :   std::vector<unsigned> shape(2);
     118         139 :   shape[0]=nvals;
     119         139 :   shape[1]=getNumberOfArguments();
     120         139 :   addValue( shape );
     121         139 :   if( periodic ) {
     122          24 :     setPeriodic( smin, smax );
     123             :   } else {
     124         115 :     setNotPeriodic();
     125             :   }
     126             :   // And store this value
     127         139 :   getPntrToComponent(0)->buildDataStore();
     128         139 :   getPntrToComponent(0)->reshapeMatrixStore( shape[1] );
     129             :   // Setup everything so we can build the store
     130         139 :   done_in_chain=true;
     131         139 :   ActionWithVector* av=dynamic_cast<ActionWithVector*>( getPntrToArgument(0)->getPntrToAction() );
     132         139 :   if( av ) {
     133          99 :     const ActionWithVector* head0 = av->getFirstActionInChain();
     134         996 :     for(unsigned i=0; i<getNumberOfArguments(); ++i) {
     135         928 :       ActionWithVector* avv=dynamic_cast<ActionWithVector*>( getPntrToArgument(i)->getPntrToAction() );
     136         928 :       if( !avv ) {
     137           4 :         continue;
     138             :       }
     139         924 :       if( head0!=avv->getFirstActionInChain() ) {
     140          31 :         done_in_chain=false;
     141          31 :         break;
     142             :       }
     143             :     }
     144             :   } else {
     145          40 :     done_in_chain=false;
     146             :   }
     147         139 :   unsigned nder = buildArgumentStore(0);
     148             :   // This checks which values have been stored
     149         139 :   stored.resize( getNumberOfArguments() );
     150         139 :   std::string headstr=getFirstActionInChain()->getLabel();
     151        1228 :   for(unsigned i=0; i<stored.size(); ++i) {
     152        1089 :     stored[i] = getPntrToArgument(i)->ignoreStoredValue( headstr );
     153             :   }
     154         139 : }
     155             : 
     156          23 : void VStack::getMatrixColumnTitles( std::vector<std::string>& argnames ) const {
     157          60 :   for(unsigned j=0; j<getNumberOfArguments(); ++j) {
     158          37 :     if( (getPntrToArgument(j)->getPntrToAction())->getName()=="COLLECT" ) {
     159          17 :       ActionWithArguments* aa = dynamic_cast<ActionWithArguments*>( getPntrToArgument(j)->getPntrToAction() );
     160          17 :       plumed_assert( aa && aa->getNumberOfArguments()==1 );
     161          17 :       argnames.push_back( (aa->getPntrToArgument(0))->getName() );
     162             :     } else {
     163          20 :       argnames.push_back( getPntrToArgument(j)->getName() );
     164             :     }
     165             :   }
     166          23 : }
     167             : 
     168        3100 : void VStack::prepare() {
     169        3100 :   ActionWithVector::prepare();
     170        3100 :   if( getPntrToArgument(0)->getRank()==0 || getPntrToArgument(0)->getShape()[0]==getPntrToComponent(0)->getShape()[0] ) {
     171        3082 :     return ;
     172             :   }
     173          18 :   std::vector<unsigned> shape(2);
     174          18 :   shape[0] = getPntrToArgument(0)->getShape()[0];
     175          18 :   shape[1] = getNumberOfArguments();
     176          18 :   getPntrToComponent(0)->setShape(shape);
     177          18 :   getPntrToComponent(0)->reshapeMatrixStore( shape[1] );
     178             : }
     179             : 
     180      214418 : void VStack::setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const {
     181      214418 :   unsigned nargs = getNumberOfArguments();
     182      214418 :   unsigned nvals = getConstPntrToComponent(0)->getShape()[0];
     183      214418 :   if( indices.size()!=nargs+1 ) {
     184       51622 :     indices.resize( nargs+1 );
     185             :   }
     186     4287565 :   for(unsigned i=0; i<nargs; ++i) {
     187     4073147 :     indices[i+1] = nvals + i;
     188             :   }
     189             :   myvals.setSplitIndex( nargs + 1 );
     190      214418 : }
     191             : 
     192     4073147 : void VStack::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
     193     4073147 :   unsigned ind2 = index2;
     194     4073147 :   if( index2>=getConstPntrToComponent(0)->getShape()[0] ) {
     195     4073147 :     ind2 = index2 - getConstPntrToComponent(0)->getShape()[0];
     196             :   }
     197     4073147 :   myvals.addValue( getConstPntrToComponent(0)->getPositionInStream(), getArgumentElement( ind2, index1, myvals ) );
     198             : 
     199     4073147 :   if( doNotCalculateDerivatives() ) {
     200      380548 :     return;
     201             :   }
     202     3692599 :   addDerivativeOnVectorArgument( stored[ind2], 0, ind2, index1, 1.0, myvals );
     203             : }
     204             : 
     205      214418 : void VStack::runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const {
     206      214418 :   if( doNotCalculateDerivatives() || !matrixChainContinues() ) {
     207             :     return ;
     208             :   }
     209             : 
     210           3 :   unsigned nmat = getConstPntrToComponent(0)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
     211             :   std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) );
     212           3 :   plumed_assert( nmat_ind<matrix_indices.size() );
     213          12 :   for(unsigned i=0; i<getNumberOfArguments(); ++i) {
     214             :     bool found=false;
     215           9 :     ActionWithValue* iav = getPntrToArgument(i)->getPntrToAction();
     216           9 :     for(unsigned j=0; j<i; ++j) {
     217           6 :       if( iav==getPntrToArgument(j)->getPntrToAction() ) {
     218             :         found=true;
     219             :         break;
     220             :       }
     221             :     }
     222           9 :     if( found ) {
     223           6 :       continue ;
     224             :     }
     225             : 
     226             :     unsigned istrn = getPntrToArgument(i)->getPositionInStream();
     227          48 :     for(unsigned k=0; k<myvals.getNumberActive(istrn); ++k) {
     228          45 :       matrix_indices[nmat_ind] = myvals.getActiveIndex(istrn,k);
     229          45 :       nmat_ind++;
     230             :     }
     231             :   }
     232             :   myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind );
     233             : }
     234             : 
     235             : }
     236             : }