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          96 : void MatrixTimesMatrix::registerKeywords( Keywords& keys ) {
      69          96 :   ActionWithMatrix::registerKeywords(keys);
      70         192 :   keys.addInputKeyword("compulsory","ARG","matrix","the label of the two matrices from which the product is calculated");
      71          96 :   keys.addFlag("SQUARED",false,"calculate the squares of the dissimilarities (this option cannot be used with MATRIX_PRODUCT)");
      72         192 :   keys.setValueDescription("matrix","the product of the two input matrices");
      73          96 : }
      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( myval->getShape()[0]==getPntrToArgument(0)->getShape()[0] && myval->getShape()[1]==getPntrToArgument(1)->getShape()[1] ) {
     116         155 :     return;
     117             :   }
     118          17 :   std::vector<unsigned> shape(2);
     119          17 :   shape[0]=getPntrToArgument(0)->getShape()[0];
     120          17 :   shape[1]=getPntrToArgument(1)->getShape()[1];
     121          17 :   myval->setShape(shape);
     122          17 :   if( myval->valueIsStored() ) {
     123          17 :     myval->reshapeMatrixStore( shape[1] );
     124             :   }
     125             : }
     126             : 
     127           0 : void MatrixTimesMatrix::getAdditionalTasksRequired( ActionWithVector* action, std::vector<unsigned>& atasks ) {
     128             : 
     129           0 :   ActionWithMatrix* adj=dynamic_cast<ActionWithMatrix*>( getPntrToArgument(0)->getPntrToAction() );
     130           0 :   if( !adj->isAdjacencyMatrix() ) {
     131             :     return;
     132             :   }
     133           0 :   adj->retrieveAtoms();
     134           0 :   adj->getAdditionalTasksRequired( action, atasks );
     135             : }
     136             : 
     137        3233 : void MatrixTimesMatrix::setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const {
     138        3233 :   unsigned start_n = getPntrToArgument(0)->getShape()[0], size_v = getPntrToArgument(1)->getShape()[1];
     139        3233 :   if( indices.size()!=size_v+1 ) {
     140         175 :     indices.resize( size_v+1 );
     141             :   }
     142      560857 :   for(unsigned i=0; i<size_v; ++i) {
     143      557624 :     indices[i+1] = start_n + i;
     144             :   }
     145             :   myvals.setSplitIndex( size_v + 1 );
     146        3233 : }
     147             : 
     148      926514 : void MatrixTimesMatrix::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
     149      926514 :   unsigned ostrn = getConstPntrToComponent(0)->getPositionInStream(), ind2=index2;
     150      926514 :   if( index2>=getPntrToArgument(0)->getShape()[0] ) {
     151      654964 :     ind2 = index2 - getPntrToArgument(0)->getShape()[0];
     152             :   }
     153             : 
     154             :   Value* myarg = getPntrToArgument(0);
     155             :   unsigned nmult=myarg->getRowLength(index1);
     156             :   double matval=0;
     157      926514 :   std::vector<double>  dvec1(nmult), dvec2(nmult);
     158    12593462 :   for(unsigned i=0; i<nmult; ++i) {
     159    11666948 :     unsigned kind = myarg->getRowIndex( index1, i );
     160    11666948 :     double val1 = getElementOfMatrixArgument( 0, index1, kind, myvals );
     161    11666948 :     double val2 = getElementOfMatrixArgument( 1, kind, ind2, myvals );
     162    11666948 :     if( getName()=="DISSIMILARITIES" ) {
     163     1637093 :       double tmp = getPntrToArgument(0)->difference(val2, val1);
     164     1637093 :       matval += tmp*tmp;
     165     1637093 :       if( !squared ) {
     166         763 :         dvec1[i] = 2*tmp;
     167         763 :         dvec2[i] = -2*tmp;
     168    10575108 :         continue;
     169             :       } else {
     170     1636330 :         val2 = -2*tmp;
     171     1636330 :         val1 = 2*tmp;
     172             :       }
     173             :     } else {
     174    10029855 :       matval+= val1*val2;
     175             :     }
     176             : 
     177    11666185 :     if( doNotCalculateDerivatives() ) {
     178    10574345 :       continue;
     179             :     }
     180             : 
     181     1091840 :     addDerivativeOnMatrixArgument( stored_matrix1, 0, 0, index1, kind, val2, myvals );
     182     1091840 :     addDerivativeOnMatrixArgument( stored_matrix2, 0, 1, kind, ind2, val1, myvals );
     183             :   }
     184             :   // And add this part of the product
     185      926514 :   if( !squared ) {
     186         388 :     matval = sqrt(matval);
     187             :   }
     188      926514 :   myvals.addValue( ostrn, matval );
     189      926514 :   if( squared || doNotCalculateDerivatives() ) {
     190             :     return;
     191             :   }
     192             : 
     193           0 :   for(unsigned i=0; i<nmult; ++i) {
     194           0 :     unsigned kind = myarg->getRowIndex( index1, i );
     195           0 :     addDerivativeOnMatrixArgument( stored_matrix1, 0, 0, index1, kind, dvec1[i]/(2*matval), myvals );
     196           0 :     addDerivativeOnMatrixArgument( stored_matrix2, 0, 1, kind, ind2, dvec2[i]/(2*matval), myvals );
     197             :   }
     198             : }
     199             : 
     200        8116 : void MatrixTimesMatrix::runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const {
     201        8116 :   if( doNotCalculateDerivatives() || !matrixChainContinues() ) {
     202             :     return ;
     203             :   }
     204             : 
     205        1734 :   unsigned mat1s = ival*getPntrToArgument(0)->getShape()[1];
     206        1734 :   unsigned nmult = getPntrToArgument(0)->getShape()[1], ss = getPntrToArgument(1)->getShape()[1];
     207        1734 :   unsigned nmat = getConstPntrToComponent(0)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
     208             :   std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) );
     209             :   unsigned ntwo_atoms = myvals.getSplitIndex();
     210       21092 :   for(unsigned j=0; j<nmult; ++j) {
     211       19358 :     matrix_indices[nmat_ind] = mat1s + j;
     212       19358 :     nmat_ind++;
     213     1098372 :     for(unsigned i=1; i<ntwo_atoms; ++i) {
     214     1079014 :       unsigned ind2 = indices[i];
     215     1079014 :       if( ind2>=getPntrToArgument(0)->getShape()[0] ) {
     216      299524 :         ind2 = indices[i] - getPntrToArgument(0)->getShape()[0];
     217             :       }
     218     1079014 :       matrix_indices[nmat_ind] = arg_deriv_starts[1] + j*ss + ind2;
     219     1079014 :       nmat_ind++;
     220             :     }
     221             :   }
     222             :   myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind );
     223             : }
     224             : 
     225             : }
     226             : }