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 "MatrixOperationBase.h"
      23             : #include "core/ActionRegister.h"
      24             : 
      25             : //+PLUMEDOC MCOLVAR TRANSPOSE
      26             : /*
      27             : Calculate the transpose of a matrix
      28             : 
      29             : This action takes a matrix in input and calculates the input matrix's [tranpose](https://en.wikipedia.org/wiki/Transpose).
      30             : The following example shows how you can use this to calculate coordination numbers of species A with species B and vice
      31             : versa.
      32             : 
      33             : ```plumed
      34             : # Calculate the contact matrix between the two groups
      35             : c1: CONTACT_MATRIX GROUPA=1-10 GROUPB=11-30 SWITCH={RATIONAL R_0=0.1}
      36             : # Calculate the cooordination numbers for the atoms in group A by multiplying by a vector of ones
      37             : onesB: ONES SIZE=20
      38             : coordA: MATRIX_VECTOR_PRODUCT ARG=c1,onesB
      39             : # Transpose the contact matrix
      40             : c1T: TRANSPOSE ARG=c1
      41             : # Calculate the coordination number for the atoms in group B by multiplying the transpose by a vector of ones
      42             : onesA: ONES SIZE=10
      43             : coordB: MATRIX_VECTOR_PRODUCT ARG=c1T,onesA
      44             : # Output the two vectors of coordination numbers to a file
      45             : PRINT ARG=coordA,coordB FILE=colvar
      46             : ```
      47             : 
      48             : Another useful example where the transpose can be used is shown below.  In this input the [DISTANCE](DISTANCE.md) command
      49             : is used to calculate the orientation of a collection of molecules.  We then can then use the [VSTACK](VSTACK.md), TRANSPOSE and the
      50             : [MATRIX_PRODUCT](MATRIX_PRODUCT.md) commands to calculate the dot products between all these vectors
      51             : 
      52             : ```plumed
      53             : # Calculate the vectors connecting these three pairs of atoms
      54             : d: DISTANCE COMPONENTS ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6
      55             : # Construct a matrix that contains all the components of the vectors calculated
      56             : v: VSTACK ARG=d.x,d.y,d.z
      57             : # Transpose v
      58             : vT: TRANSPOSE ARG=v
      59             : # And now calculate the 3x3 matrix of dot products
      60             : m: MATRIX_PRODUCT ARG=v,vT
      61             : # And output the matrix product to a file
      62             : PRINT ARG=m FILE=colvar
      63             : ```
      64             : 
      65             : */
      66             : //+ENDPLUMEDOC
      67             : 
      68             : namespace PLMD {
      69             : namespace matrixtools {
      70             : 
      71             : class TransposeMatrix : public MatrixOperationBase {
      72             : public:
      73             :   static void registerKeywords( Keywords& keys );
      74             : /// Constructor
      75             :   explicit TransposeMatrix(const ActionOptions&);
      76             : ///
      77         256 :   unsigned getNumberOfDerivatives() override {
      78         256 :     return 0;
      79             :   }
      80             : ///
      81             :   void prepare() override ;
      82             : ///
      83             :   void calculate() override ;
      84             : ///
      85             :   void apply() override ;
      86             : ///
      87             :   double getForceOnMatrixElement( const unsigned& jrow, const unsigned& krow ) const override;
      88             : };
      89             : 
      90             : PLUMED_REGISTER_ACTION(TransposeMatrix,"TRANSPOSE")
      91             : 
      92         279 : void TransposeMatrix::registerKeywords( Keywords& keys ) {
      93         279 :   MatrixOperationBase::registerKeywords( keys );
      94         558 :   keys.addInputKeyword("compulsory","ARG","vector/matrix","the label of the vector or matrix that should be transposed");
      95         558 :   keys.setValueDescription("vector/matrix","the transpose of the input matrix");
      96         279 : }
      97             : 
      98         162 : TransposeMatrix::TransposeMatrix(const ActionOptions& ao):
      99             :   Action(ao),
     100         162 :   MatrixOperationBase(ao) {
     101         162 :   if( getPntrToArgument(0)->isSymmetric() ) {
     102           0 :     error("input matrix is symmetric.  Transposing will achieve nothing!");
     103             :   }
     104             :   std::vector<unsigned> shape;
     105         162 :   if( getPntrToArgument(0)->getRank()==0 ) {
     106           0 :     error("transposing a scalar?");
     107         162 :   } else if( getPntrToArgument(0)->getRank()==1 ) {
     108          17 :     shape.resize(2);
     109          17 :     shape[0]=1;
     110          17 :     shape[1]=getPntrToArgument(0)->getShape()[0];
     111         145 :   } else if( getPntrToArgument(0)->getShape()[0]==1 ) {
     112          61 :     shape.resize(1);
     113          61 :     shape[0] = getPntrToArgument(0)->getShape()[1];
     114             :   } else {
     115          84 :     shape.resize(2);
     116          84 :     shape[0]=getPntrToArgument(0)->getShape()[1];
     117          84 :     shape[1]=getPntrToArgument(0)->getShape()[0];
     118             :   }
     119         162 :   addValue( shape );
     120         162 :   if( getPntrToArgument(0)->isPeriodic() ) {
     121             :     std::string smin, smax;
     122          24 :     getPntrToArgument(0)->getDomain( smin, smax );
     123          24 :     setPeriodic( smin, smax );
     124             :   } else {
     125         138 :     setNotPeriodic();
     126             :   }
     127         162 :   getPntrToComponent(0)->buildDataStore();
     128         162 :   if( shape.size()==2 ) {
     129         101 :     getPntrToComponent(0)->reshapeMatrixStore( shape[1] );
     130             :   }
     131         162 : }
     132             : 
     133        4821 : void TransposeMatrix::prepare() {
     134        4821 :   Value* myval = getPntrToComponent(0);
     135             :   Value* myarg = getPntrToArgument(0);
     136        4821 :   if( myarg->getRank()==1 ) {
     137         586 :     if( myval->getShape()[0]!=1 || myval->getShape()[1]!=myarg->getShape()[0] ) {
     138           6 :       std::vector<unsigned> shape(2);
     139           6 :       shape[0] = 1;
     140           6 :       shape[1] = myarg->getShape()[0];
     141           6 :       myval->setShape( shape );
     142           6 :       myval->reshapeMatrixStore( shape[1] );
     143             :     }
     144        4235 :   } else if( myarg->getShape()[0]==1 ) {
     145        2392 :     if( myval->getShape()[0]!=myarg->getShape()[1] ) {
     146           6 :       std::vector<unsigned> shape(1);
     147           6 :       shape[0] = myarg->getShape()[1];
     148           6 :       myval->setShape( shape );
     149             :     }
     150        1843 :   } else if( myarg->getShape()[0]!=myval->getShape()[1] || myarg->getShape()[1]!=myval->getShape()[0] ) {
     151          19 :     std::vector<unsigned> shape(2);
     152          19 :     shape[0] = myarg->getShape()[1];
     153          19 :     shape[1] = myarg->getShape()[0];
     154          19 :     myval->setShape( shape );
     155          19 :     myval->reshapeMatrixStore( shape[1] );
     156             :   }
     157        4821 : }
     158             : 
     159        4225 : void TransposeMatrix::calculate() {
     160             :   // Retrieve the non-zero pairs
     161             :   Value* myarg=getPntrToArgument(0);
     162        4225 :   Value* myval=getPntrToComponent(0);
     163        4225 :   if( myarg->getRank()<=1 || myval->getRank()==1 ) {
     164        2389 :     if( myarg->getRank()<=1 && myval->getShape()[1]!=myarg->getShape()[0] ) {
     165           0 :       std::vector<unsigned> shape( 2 );
     166           0 :       shape[0] = 1;
     167           0 :       shape[1] = myarg->getShape()[0];
     168           0 :       myval->setShape( shape );
     169           0 :       myval->reshapeMatrixStore( shape[1] );
     170        2389 :     } else if( myval->getRank()==1 && myval->getShape()[0]!=myarg->getShape()[1] ) {
     171           0 :       std::vector<unsigned> shape( 1 );
     172           0 :       shape[0] = myarg->getShape()[1];
     173           0 :       myval->setShape( shape );
     174             :     }
     175        2389 :     unsigned nv=myarg->getNumberOfValues();
     176       49015 :     for(unsigned i=0; i<nv; ++i) {
     177       46626 :       myval->set( i, myarg->get(i) );
     178             :     }
     179             :   } else {
     180        1836 :     if( myarg->getShape()[0]!=myval->getShape()[1] || myarg->getShape()[1]!=myval->getShape()[0] ) {
     181           0 :       std::vector<unsigned> shape( 2 );
     182           0 :       shape[0] = myarg->getShape()[1];
     183           0 :       shape[1] = myarg->getShape()[0];
     184           0 :       myval->setShape( shape );
     185           0 :       myval->reshapeMatrixStore( shape[1] );
     186             :     }
     187             :     std::vector<double> vals;
     188             :     std::vector<std::pair<unsigned,unsigned> > pairs;
     189        1836 :     std::vector<unsigned> shape( myval->getShape() );
     190        1836 :     unsigned nedge=0;
     191        1836 :     myarg->retrieveEdgeList( nedge, pairs, vals );
     192     2758860 :     for(unsigned i=0; i<nedge; ++i) {
     193     2757024 :       myval->set( pairs[i].second*shape[1] + pairs[i].first, vals[i] );
     194             :     }
     195             :   }
     196        4225 : }
     197             : 
     198        4144 : void TransposeMatrix::apply() {
     199        4144 :   if( doNotCalculateDerivatives() ) {
     200             :     return;
     201             :   }
     202             : 
     203             :   // Apply force on the matrix
     204        1930 :   if( getPntrToComponent(0)->forcesWereAdded() ) {
     205             :     Value* myarg=getPntrToArgument(0);
     206        1930 :     Value* myval=getPntrToComponent(0);
     207        1930 :     if( myarg->getRank()<=1 || myval->getRank()==1 ) {
     208         588 :       unsigned nv=myarg->getNumberOfValues();
     209        2408 :       for(unsigned i=0; i<nv; ++i) {
     210        1820 :         myarg->addForce( i, myval->getForce(i) );
     211             :       }
     212             :     } else {
     213        1342 :       MatrixOperationBase::apply();
     214             :     }
     215             :   }
     216             : }
     217             : 
     218     4124572 : double TransposeMatrix::getForceOnMatrixElement( const unsigned& jrow, const unsigned& kcol ) const {
     219     4124572 :   return getConstPntrToComponent(0)->getForce(kcol*getConstPntrToComponent(0)->getShape()[1]+jrow);
     220             : }
     221             : 
     222             : 
     223             : }
     224             : }