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             : private:
      73             : /// The array for matrix bookeeping
      74             :   std::vector<unsigned> matrix_bookeeping;
      75             : /// The locations of each element of the transposed matrix in the original matrix
      76             :   std::vector<unsigned> force_mapper;
      77             : public:
      78             :   static void registerKeywords( Keywords& keys );
      79             : /// Constructor
      80             :   explicit TransposeMatrix(const ActionOptions&);
      81             : ///
      82         284 :   unsigned getNumberOfDerivatives() override {
      83         284 :     return 0;
      84             :   }
      85             : ///
      86             :   void prepare() override ;
      87             : ///
      88             :   void calculate() override ;
      89             : ///
      90             :   void apply() override ;
      91             : ///
      92           0 :   double getForceOnMatrixElement( const unsigned& jrow, const unsigned& krow ) const override {
      93           0 :     plumed_error();
      94             :   }
      95             : };
      96             : 
      97             : PLUMED_REGISTER_ACTION(TransposeMatrix,"TRANSPOSE")
      98             : 
      99         319 : void TransposeMatrix::registerKeywords( Keywords& keys ) {
     100         319 :   MatrixOperationBase::registerKeywords( keys );
     101         638 :   keys.addInputKeyword("compulsory","ARG","vector/matrix","the label of the vector or matrix that should be transposed");
     102         638 :   keys.setValueDescription("vector/matrix","the transpose of the input matrix");
     103         319 :   keys.add("hidden","MASKED_INPUT_ALLOWED","turns on that you are allowed to use masked inputs ");
     104         319 : }
     105             : 
     106         184 : TransposeMatrix::TransposeMatrix(const ActionOptions& ao):
     107             :   Action(ao),
     108         184 :   MatrixOperationBase(ao) {
     109         184 :   if( getPntrToArgument(0)->isSymmetric() ) {
     110           0 :     error("input matrix is symmetric.  Transposing will achieve nothing!");
     111             :   }
     112         184 :   std::vector<std::size_t> shape;
     113         184 :   if( getPntrToArgument(0)->getRank()==0 ) {
     114           0 :     error("transposing a scalar?");
     115         184 :   } else if( getPntrToArgument(0)->getRank()==1 ) {
     116          17 :     shape.resize(2);
     117          17 :     shape[0]=1;
     118          17 :     shape[1]=getPntrToArgument(0)->getShape()[0];
     119         167 :   } else if( getPntrToArgument(0)->getShape()[0]==1 ) {
     120          61 :     shape.resize(1);
     121          61 :     shape[0] = getPntrToArgument(0)->getShape()[1];
     122             :   } else {
     123         106 :     shape.resize(2);
     124         106 :     shape[0]=getPntrToArgument(0)->getShape()[1];
     125         106 :     shape[1]=getPntrToArgument(0)->getShape()[0];
     126             :   }
     127         184 :   addValue( shape );
     128         184 :   if( getPntrToArgument(0)->isPeriodic() ) {
     129             :     std::string smin, smax;
     130          24 :     getPntrToArgument(0)->getDomain( smin, smax );
     131          24 :     setPeriodic( smin, smax );
     132             :   } else {
     133         160 :     setNotPeriodic();
     134             :   }
     135         184 :   if( getPntrToArgument(0)->isDerivativeZeroWhenValueIsZero() ) {
     136          68 :     getPntrToComponent(0)->setDerivativeIsZeroWhenValueIsZero();
     137             :   }
     138         184 : }
     139             : 
     140        5518 : void TransposeMatrix::prepare() {
     141        5518 :   Value* myval = getPntrToComponent(0);
     142             :   Value* myarg = getPntrToArgument(0);
     143        5518 :   if( myarg->getRank()==1 ) {
     144         586 :     if( myval->getShape()[0]!=1 || myval->getShape()[1]!=myarg->getShape()[0] ) {
     145           6 :       std::vector<std::size_t> shape(2);
     146           6 :       shape[0] = 1;
     147           6 :       shape[1] = myarg->getShape()[0];
     148           6 :       myval->setShape( shape );
     149           6 :     }
     150         586 :     if( myval->getNumberOfColumns()!=myarg->getShape()[0] ) {
     151          13 :       myval->reshapeMatrixStore( myarg->getShape()[0]  );
     152             :     }
     153        4932 :   } else if( myarg->getShape()[0]==1 ) {
     154        2392 :     if( myval->getShape()[0]!=myarg->getShape()[1] ) {
     155           6 :       std::vector<std::size_t> shape(1);
     156           6 :       shape[0] = myarg->getShape()[1];
     157           6 :       myval->setShape( shape );
     158           6 :     }
     159        2540 :   } else if( myarg->getShape()[0]!=myval->getShape()[1] || myarg->getShape()[1]!=myval->getShape()[0] ) {
     160          19 :     std::vector<std::size_t> shape(2);
     161          19 :     shape[0] = myarg->getShape()[1];
     162          19 :     shape[1] = myarg->getShape()[0];
     163          19 :     myval->setShape( shape );
     164          19 :   }
     165        5518 : }
     166             : 
     167        4924 : void TransposeMatrix::calculate() {
     168             :   // Retrieve the non-zero pairs
     169             :   Value* myarg=getPntrToArgument(0);
     170        4924 :   Value* myval=getPntrToComponent(0);
     171        4924 :   if( myarg->getRank()<=1 || myval->getRank()==1 ) {
     172        2389 :     unsigned nv=myarg->getNumberOfValues();
     173       49015 :     for(unsigned i=0; i<nv; ++i) {
     174       46626 :       myval->set( i, myarg->get(i) );
     175             :     }
     176             :   } else {
     177             :     // Find the longest column
     178        2535 :     std::size_t maxcol = myarg->getLengthOfLongestColumn();
     179        2535 :     myval->reshapeMatrixStore( maxcol );
     180             :     // the store keeps at least one column (see Value::reshapeMatrixStore)
     181             :     maxcol = myval->getNumberOfColumns();
     182        2535 :     matrix_bookeeping.resize( myval->getShape()[0]*(1+maxcol) );
     183        2535 :     force_mapper.resize( myval->getShape()[0]*maxcol );
     184      141089 :     for(unsigned i=0; i<myval->getShape()[0]; ++i) {
     185      138554 :       matrix_bookeeping[i*(1+maxcol)] = 0;
     186             :     }
     187             :     // Now get the bookeeping
     188        2535 :     unsigned arg_ncol = myarg->getNumberOfColumns();
     189      330784 :     for(unsigned i=0; i<myarg->getShape()[0]; ++i) {
     190             :       unsigned nr = myarg->getRowLength(i);
     191     6611398 :       for(unsigned j=0; j<nr; ++j) {
     192     6283149 :         unsigned ind = myarg->getRowIndex( i, j );
     193     6283149 :         unsigned startpos = ind*(1+maxcol);
     194     6283149 :         matrix_bookeeping[startpos + 1 + matrix_bookeeping[startpos]] = i;
     195     6283149 :         myval->set( ind*maxcol + matrix_bookeeping[startpos], myarg->get( i*arg_ncol+j, false)  );
     196     6283149 :         force_mapper[ind*maxcol + matrix_bookeeping[startpos]] = i*arg_ncol+j;
     197     6283149 :         matrix_bookeeping[startpos]++;
     198             :       }
     199             :     }
     200             :     // And copy all the bookeeping to the output value
     201     6839218 :     for(unsigned i=0; i<matrix_bookeeping.size(); ++i) {
     202             :       myval->setMatrixBookeepingElement( i, matrix_bookeeping[i] );
     203             :     }
     204             :   }
     205        4924 : }
     206             : 
     207        4843 : void TransposeMatrix::apply() {
     208        4843 :   if( doNotCalculateDerivatives() ) {
     209             :     return;
     210             :   }
     211             : 
     212             :   // Apply force on the matrix
     213        2572 :   if( getPntrToComponent(0)->forcesWereAdded() ) {
     214             :     Value* myarg=getPntrToArgument(0);
     215        2551 :     Value* myval=getPntrToComponent(0);
     216        2551 :     if( myarg->getRank()<=1 || myval->getRank()==1 ) {
     217         588 :       unsigned nv=myarg->getNumberOfValues();
     218        2408 :       for(unsigned i=0; i<nv; ++i) {
     219        1820 :         myarg->addForce( i, myval->getForce(i) );
     220             :       }
     221             :     } else {
     222        1963 :       unsigned nval_cols = myval->getNumberOfColumns();
     223       49170 :       for(unsigned i=0; i<myval->getShape()[0]; ++i) {
     224             :         unsigned nr = myval->getRowLength(i);
     225     1935958 :         for(unsigned j=0; j<nr; ++j) {
     226     1888751 :           myarg->addForce( force_mapper[nval_cols*i+j], myval->getForce(nval_cols*i + j), false );
     227             :         }
     228             :       }
     229             :     }
     230             :   }
     231             : }
     232             : 
     233             : }
     234             : }