Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2011-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 "core/ActionWithMatrix.h"
      23             : #include "core/ActionRegister.h"
      24             : 
      25             : //+PLUMEDOC MCOLVAR MATRIX_PRODUCT
      26             : /*
      27             : Calculate the product of two matrices
      28             : 
      29             : \par Examples
      30             : 
      31             : */
      32             : //+ENDPLUMEDOC
      33             : 
      34             : //+PLUMEDOC ANALYSIS DISSIMILARITIES
      35             : /*
      36             : Calculate the matrix of dissimilarities between a trajectory of atomic configurations.
      37             : 
      38             : \par Examples
      39             : 
      40             : */
      41             : //+ENDPLUMEDOC
      42             : 
      43             : namespace PLMD {
      44             : namespace matrixtools {
      45             : 
      46             : class MatrixTimesMatrix : public ActionWithMatrix {
      47             : private:
      48             :   bool squared;
      49             :   unsigned nderivatives;
      50             :   bool stored_matrix1, stored_matrix2;
      51             : public:
      52             :   static void registerKeywords( Keywords& keys );
      53             :   explicit MatrixTimesMatrix(const ActionOptions&);
      54             :   void prepare() override ;
      55             :   unsigned getNumberOfDerivatives();
      56        1297 :   unsigned getNumberOfColumns() const override {
      57        1297 :     return getConstPntrToComponent(0)->getShape()[1];
      58             :   }
      59             :   void getAdditionalTasksRequired( ActionWithVector* action, std::vector<unsigned>& atasks ) override ;
      60             :   void setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const ;
      61             :   void performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const override;
      62             :   void runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const override ;
      63             : };
      64             : 
      65             : PLUMED_REGISTER_ACTION(MatrixTimesMatrix,"MATRIX_PRODUCT")
      66             : PLUMED_REGISTER_ACTION(MatrixTimesMatrix,"DISSIMILARITIES")
      67             : 
      68         100 : void MatrixTimesMatrix::registerKeywords( Keywords& keys ) {
      69         100 :   ActionWithMatrix::registerKeywords(keys);
      70         100 :   keys.use("ARG");
      71         200 :   keys.addFlag("SQUARED",false,"calculate the squares of the dissimilarities (this option cannot be used with MATRIX_PRODUCT)");
      72         100 :   keys.setValueDescription("the product of the two input matrices");
      73         100 : }
      74             : 
      75          57 : MatrixTimesMatrix::MatrixTimesMatrix(const ActionOptions&ao):
      76             :   Action(ao),
      77          57 :   ActionWithMatrix(ao) {
      78          57 :   if( getNumberOfArguments()!=2 ) {
      79           0 :     error("should be two arguments to this action, a matrix and a vector");
      80             :   }
      81          57 :   if( getPntrToArgument(0)->getRank()!=2 || getPntrToArgument(0)->hasDerivatives() ) {
      82           0 :     error("first argument to this action should be a matrix");
      83             :   }
      84          57 :   if( getPntrToArgument(1)->getRank()!=2 || getPntrToArgument(1)->hasDerivatives() ) {
      85           0 :     error("second argument to this action should be a matrix");
      86             :   }
      87          57 :   if( getPntrToArgument(0)->getShape()[1]!=getPntrToArgument(1)->getShape()[0] ) {
      88           0 :     error("number of columns in first matrix does not equal number of rows in second matrix");
      89             :   }
      90          57 :   std::vector<unsigned> shape(2);
      91          57 :   shape[0]=getPntrToArgument(0)->getShape()[0];
      92          57 :   shape[1]=getPntrToArgument(1)->getShape()[1];
      93          57 :   addValue( shape );
      94          57 :   setNotPeriodic();
      95          57 :   nderivatives = buildArgumentStore(0);
      96          57 :   std::string headstr=getFirstActionInChain()->getLabel();
      97          57 :   stored_matrix1 = getPntrToArgument(0)->ignoreStoredValue( headstr );
      98          57 :   stored_matrix2 = getPntrToArgument(1)->ignoreStoredValue( headstr );
      99          57 :   if( getName()=="DISSIMILARITIES" ) {
     100          12 :     parseFlag("SQUARED",squared);
     101          12 :     if( squared ) {
     102           7 :       log.printf("  calculating the squares of the dissimilarities \n");
     103             :     }
     104             :   } else {
     105          45 :     squared=true;
     106             :   }
     107          57 : }
     108             : 
     109          65 : unsigned MatrixTimesMatrix::getNumberOfDerivatives() {
     110          65 :   return nderivatives;
     111             : }
     112             : 
     113         172 : void MatrixTimesMatrix::prepare() {
     114         172 :   Value* myval = getPntrToComponent(0);
     115         172 :   if( !getPntrToComponent(0)->isDerivativeZeroWhenValueIsZero() ) {
     116         172 :     if( getPntrToArgument(0)->getNumberOfColumns()>0 && getPntrToArgument(0)->getNumberOfColumns()<getPntrToArgument(0)->getShape()[1] ) {
     117           0 :       if( !doNotCalculateDerivatives() ) {
     118           0 :         error("cannot calculate derivatives for this action with sparse matrices");
     119           0 :       } else if( getName()=="DISSIMILARITIES" ) {
     120           0 :         error("cannot calculate dissimilarities for sparse matrices");
     121             :       }
     122             :     }
     123         172 :     if(getPntrToArgument(1)->getNumberOfColumns()>0 && getPntrToArgument(1)->getNumberOfColumns()<getPntrToArgument(1)->getShape()[1] ) {
     124           0 :       if( !doNotCalculateDerivatives() ) {
     125           0 :         error("cannot calculate derivatives for this action with sparse matrices");
     126           0 :       } else if( getName()=="DISSIMILARITIES" ) {
     127           0 :         error("cannot calculate dissimilarities for sparse matrices");
     128             :       }
     129             :     }
     130             :   }
     131         172 :   if( myval->getShape()[0]==getPntrToArgument(0)->getShape()[0] && myval->getShape()[1]==getPntrToArgument(1)->getShape()[1] ) {
     132         155 :     return;
     133             :   }
     134          17 :   std::vector<unsigned> shape(2);
     135          17 :   shape[0]=getPntrToArgument(0)->getShape()[0];
     136          17 :   shape[1]=getPntrToArgument(1)->getShape()[1];
     137          17 :   myval->setShape(shape);
     138          17 :   if( myval->valueIsStored() ) {
     139          17 :     myval->reshapeMatrixStore( shape[1] );
     140             :   }
     141             : }
     142             : 
     143           0 : void MatrixTimesMatrix::getAdditionalTasksRequired( ActionWithVector* action, std::vector<unsigned>& atasks ) {
     144             : 
     145           0 :   ActionWithMatrix* adj=dynamic_cast<ActionWithMatrix*>( getPntrToArgument(0)->getPntrToAction() );
     146           0 :   if( !adj->isAdjacencyMatrix() ) {
     147             :     return;
     148             :   }
     149           0 :   adj->retrieveAtoms();
     150           0 :   adj->getAdditionalTasksRequired( action, atasks );
     151             : }
     152             : 
     153        3233 : void MatrixTimesMatrix::setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const {
     154        3233 :   unsigned start_n = getPntrToArgument(0)->getShape()[0], size_v = getPntrToArgument(1)->getShape()[1];
     155        3233 :   if( indices.size()!=size_v+1 ) {
     156         175 :     indices.resize( size_v+1 );
     157             :   }
     158      560857 :   for(unsigned i=0; i<size_v; ++i) {
     159      557624 :     indices[i+1] = start_n + i;
     160             :   }
     161             :   myvals.setSplitIndex( size_v + 1 );
     162        3233 : }
     163             : 
     164     4365496 : void MatrixTimesMatrix::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
     165     4365496 :   unsigned ostrn = getConstPntrToComponent(0)->getPositionInStream(), ind2=index2;
     166     4365496 :   if( index2>=getPntrToArgument(0)->getShape()[0] ) {
     167      654964 :     ind2 = index2 - getPntrToArgument(0)->getShape()[0];
     168             :   }
     169             : 
     170             :   Value* myarg = getPntrToArgument(0);
     171             :   unsigned nmult=myarg->getRowLength(index1);
     172             :   double matval=0;
     173     4365496 :   std::vector<double>  dvec1(nmult), dvec2(nmult);
     174   105365336 :   for(unsigned i=0; i<nmult; ++i) {
     175   100999840 :     unsigned kind = myarg->getRowIndex( index1, i );
     176   100999840 :     double val1 = getElementOfMatrixArgument( 0, index1, kind, myvals );
     177   100999840 :     double val2 = getElementOfMatrixArgument( 1, kind, ind2, myvals );
     178   100999840 :     if( getName()=="DISSIMILARITIES" ) {
     179     1637093 :       double tmp = getPntrToArgument(0)->difference(val2, val1);
     180     1637093 :       matval += tmp*tmp;
     181     1637093 :       if( !squared ) {
     182         763 :         dvec1[i] = 2*tmp;
     183         763 :         dvec2[i] = -2*tmp;
     184    10575108 :         continue;
     185             :       } else {
     186     1636330 :         val2 = -2*tmp;
     187     1636330 :         val1 = 2*tmp;
     188             :       }
     189             :     } else {
     190    99362747 :       matval+= val1*val2;
     191             :     }
     192             : 
     193   100999077 :     if( doNotCalculateDerivatives() ) {
     194    10574345 :       continue;
     195             :     }
     196             : 
     197    90424732 :     addDerivativeOnMatrixArgument( stored_matrix1, 0, 0, index1, kind, val2, myvals );
     198    90424732 :     addDerivativeOnMatrixArgument( stored_matrix2, 0, 1, kind, ind2, val1, myvals );
     199             :   }
     200             :   // And add this part of the product
     201     4365496 :   if( !squared ) {
     202         388 :     matval = sqrt(matval);
     203             :   }
     204     4365496 :   myvals.addValue( ostrn, matval );
     205     4365496 :   if( squared || doNotCalculateDerivatives() ) {
     206             :     return;
     207             :   }
     208             : 
     209           0 :   for(unsigned i=0; i<nmult; ++i) {
     210           0 :     unsigned kind = myarg->getRowIndex( index1, i );
     211           0 :     addDerivativeOnMatrixArgument( stored_matrix1, 0, 0, index1, kind, dvec1[i]/(2*matval), myvals );
     212           0 :     addDerivativeOnMatrixArgument( stored_matrix2, 0, 1, kind, ind2, dvec2[i]/(2*matval), myvals );
     213             :   }
     214             : }
     215             : 
     216       94805 : void MatrixTimesMatrix::runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const {
     217       94805 :   if( doNotCalculateDerivatives() || !matrixChainContinues() ) {
     218             :     return ;
     219             :   }
     220             : 
     221       88423 :   unsigned mat1s = ival*getPntrToArgument(0)->getShape()[1];
     222       88423 :   unsigned nmult = getPntrToArgument(0)->getShape()[1], ss = getPntrToArgument(1)->getShape()[1];
     223       88423 :   unsigned nmat = getConstPntrToComponent(0)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
     224             :   std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) );
     225             :   unsigned ntwo_atoms = myvals.getSplitIndex();
     226     2360415 :   for(unsigned j=0; j<nmult; ++j) {
     227     2271992 :     matrix_indices[nmat_ind] = mat1s + j;
     228     2271992 :     nmat_ind++;
     229    92683898 :     for(unsigned i=1; i<ntwo_atoms; ++i) {
     230    90411906 :       unsigned ind2 = indices[i];
     231    90411906 :       if( ind2>=getPntrToArgument(0)->getShape()[0] ) {
     232      299524 :         ind2 = indices[i] - getPntrToArgument(0)->getShape()[0];
     233             :       }
     234    90411906 :       matrix_indices[nmat_ind] = arg_deriv_starts[1] + j*ss + ind2;
     235    90411906 :       nmat_ind++;
     236             :     }
     237             :   }
     238             :   myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind );
     239             : }
     240             : 
     241             : }
     242             : }