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 { return getConstPntrToComponent(0)->getShape()[1]; }
      57             :   void getAdditionalTasksRequired( ActionWithVector* action, std::vector<unsigned>& atasks ) override ;
      58             :   void setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const ;
      59             :   void performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const override;
      60             :   void runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const override ;
      61             : };
      62             : 
      63             : PLUMED_REGISTER_ACTION(MatrixTimesMatrix,"MATRIX_PRODUCT")
      64             : PLUMED_REGISTER_ACTION(MatrixTimesMatrix,"DISSIMILARITIES")
      65             : 
      66          96 : void MatrixTimesMatrix::registerKeywords( Keywords& keys ) {
      67          96 :   ActionWithMatrix::registerKeywords(keys); keys.use("ARG");
      68         192 :   keys.addFlag("SQUARED",false,"calculate the squares of the dissimilarities (this option cannot be used with MATRIX_PRODUCT)");
      69          96 :   keys.setValueDescription("the product of the two input matrices");
      70          96 : }
      71             : 
      72          57 : MatrixTimesMatrix::MatrixTimesMatrix(const ActionOptions&ao):
      73             :   Action(ao),
      74          57 :   ActionWithMatrix(ao)
      75             : {
      76          57 :   if( getNumberOfArguments()!=2 ) error("should be two arguments to this action, a matrix and a vector");
      77          57 :   if( getPntrToArgument(0)->getRank()!=2 || getPntrToArgument(0)->hasDerivatives() ) error("first argument to this action should be a matrix");
      78          57 :   if( getPntrToArgument(1)->getRank()!=2 || getPntrToArgument(1)->hasDerivatives() ) error("second argument to this action should be a matrix");
      79          57 :   if( getPntrToArgument(0)->getShape()[1]!=getPntrToArgument(1)->getShape()[0] ) error("number of columns in first matrix does not equal number of rows in second matrix");
      80          57 :   std::vector<unsigned> shape(2); shape[0]=getPntrToArgument(0)->getShape()[0]; shape[1]=getPntrToArgument(1)->getShape()[1];
      81          57 :   addValue( shape ); setNotPeriodic(); nderivatives = buildArgumentStore(0);
      82          57 :   std::string headstr=getFirstActionInChain()->getLabel();
      83          57 :   stored_matrix1 = getPntrToArgument(0)->ignoreStoredValue( headstr );
      84          57 :   stored_matrix2 = getPntrToArgument(1)->ignoreStoredValue( headstr );
      85          57 :   if( getName()=="DISSIMILARITIES" ) {
      86          12 :     parseFlag("SQUARED",squared);
      87          12 :     if( squared ) log.printf("  calculating the squares of the dissimilarities \n");
      88          45 :   } else squared=true;
      89          57 : }
      90             : 
      91          65 : unsigned MatrixTimesMatrix::getNumberOfDerivatives() {
      92          65 :   return nderivatives;
      93             : }
      94             : 
      95         172 : void MatrixTimesMatrix::prepare() {
      96         172 :   Value* myval = getPntrToComponent(0);
      97         172 :   if( myval->getShape()[0]==getPntrToArgument(0)->getShape()[0] && myval->getShape()[1]==getPntrToArgument(1)->getShape()[1] ) return;
      98          17 :   std::vector<unsigned> shape(2); shape[0]=getPntrToArgument(0)->getShape()[0]; shape[1]=getPntrToArgument(1)->getShape()[1];
      99          17 :   myval->setShape(shape); if( myval->valueIsStored() ) myval->reshapeMatrixStore( shape[1] );
     100             : }
     101             : 
     102           0 : void MatrixTimesMatrix::getAdditionalTasksRequired( ActionWithVector* action, std::vector<unsigned>& atasks ) {
     103             : 
     104           0 :   ActionWithMatrix* adj=dynamic_cast<ActionWithMatrix*>( getPntrToArgument(0)->getPntrToAction() );
     105           0 :   if( !adj->isAdjacencyMatrix() ) return;
     106           0 :   adj->retrieveAtoms(); adj->getAdditionalTasksRequired( action, atasks );
     107             : }
     108             : 
     109        3233 : void MatrixTimesMatrix::setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const {
     110        3233 :   unsigned start_n = getPntrToArgument(0)->getShape()[0], size_v = getPntrToArgument(1)->getShape()[1];
     111        3233 :   if( indices.size()!=size_v+1 ) indices.resize( size_v+1 );
     112      560857 :   for(unsigned i=0; i<size_v; ++i) indices[i+1] = start_n + i;
     113             :   myvals.setSplitIndex( size_v + 1 );
     114        3233 : }
     115             : 
     116      926514 : void MatrixTimesMatrix::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
     117      926514 :   unsigned ostrn = getConstPntrToComponent(0)->getPositionInStream(), ind2=index2;
     118      926514 :   if( index2>=getPntrToArgument(0)->getShape()[0] ) ind2 = index2 - getPntrToArgument(0)->getShape()[0];
     119             : 
     120             :   Value* myarg = getPntrToArgument(0);
     121             :   unsigned nmult=myarg->getRowLength(index1); double matval=0;
     122      926514 :   std::vector<double>  dvec1(nmult), dvec2(nmult);
     123    12593462 :   for(unsigned i=0; i<nmult; ++i) {
     124    11666948 :     unsigned kind = myarg->getRowIndex( index1, i );
     125    11666948 :     double val1 = getElementOfMatrixArgument( 0, index1, kind, myvals );
     126    11666948 :     double val2 = getElementOfMatrixArgument( 1, kind, ind2, myvals );
     127    11666948 :     if( getName()=="DISSIMILARITIES" ) {
     128     1637093 :       double tmp = getPntrToArgument(0)->difference(val2, val1); matval += tmp*tmp;
     129     1637093 :       if( !squared ) {
     130    10575108 :         dvec1[i] = 2*tmp; dvec2[i] = -2*tmp; continue;
     131     1636330 :       } else { val2 = -2*tmp; val1 = 2*tmp; }
     132    10029855 :     } else matval+= val1*val2;
     133             : 
     134    11666185 :     if( doNotCalculateDerivatives() ) continue;
     135             : 
     136     1091840 :     addDerivativeOnMatrixArgument( stored_matrix1, 0, 0, index1, kind, val2, myvals );
     137     1091840 :     addDerivativeOnMatrixArgument( stored_matrix2, 0, 1, kind, ind2, val1, myvals );
     138             :   }
     139             :   // And add this part of the product
     140      926514 :   if( !squared ) matval = sqrt(matval);
     141      926514 :   myvals.addValue( ostrn, matval );
     142      926514 :   if( squared || doNotCalculateDerivatives() ) return;
     143             : 
     144           0 :   for(unsigned i=0; i<nmult; ++i) {
     145           0 :     unsigned kind = myarg->getRowIndex( index1, i );
     146           0 :     addDerivativeOnMatrixArgument( stored_matrix1, 0, 0, index1, kind, dvec1[i]/(2*matval), myvals );
     147           0 :     addDerivativeOnMatrixArgument( stored_matrix2, 0, 1, kind, ind2, dvec2[i]/(2*matval), myvals );
     148             :   }
     149             : }
     150             : 
     151        8116 : void MatrixTimesMatrix::runEndOfRowJobs( const unsigned& ival, const std::vector<unsigned> & indices, MultiValue& myvals ) const {
     152        8116 :   if( doNotCalculateDerivatives() || !matrixChainContinues() ) return ;
     153             : 
     154        1734 :   unsigned mat1s = ival*getPntrToArgument(0)->getShape()[1];
     155        1734 :   unsigned nmult = getPntrToArgument(0)->getShape()[1], ss = getPntrToArgument(1)->getShape()[1];
     156        1734 :   unsigned nmat = getConstPntrToComponent(0)->getPositionInMatrixStash(), nmat_ind = myvals.getNumberOfMatrixRowDerivatives( nmat );
     157             :   std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices( nmat ) ); unsigned ntwo_atoms = myvals.getSplitIndex();
     158       21092 :   for(unsigned j=0; j<nmult; ++j) {
     159       19358 :     matrix_indices[nmat_ind] = mat1s + j; nmat_ind++;
     160     1098372 :     for(unsigned i=1; i<ntwo_atoms; ++i) {
     161     1079014 :       unsigned ind2 = indices[i]; if( ind2>=getPntrToArgument(0)->getShape()[0] ) ind2 = indices[i] - getPntrToArgument(0)->getShape()[0];
     162     1079014 :       matrix_indices[nmat_ind] = arg_deriv_starts[1] + j*ss + ind2; nmat_ind++;
     163             :     }
     164             :   }
     165             :   myvals.setNumberOfMatrixRowDerivatives( nmat, nmat_ind );
     166             : }
     167             : 
     168             : }
     169             : }