Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 : Copyright (c) 2017-2019 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/ActionRegister.h" 23 : #include "core/ActionShortcut.h" 24 : 25 : //+PLUMEDOC REWEIGHTING COVARIANCE_MATRIX 26 : /* 27 : Calculate a covariance matix 28 : 29 : \par Examples 30 : 31 : */ 32 : //+ENDPLUMEDOC 33 : 34 : namespace PLMD { 35 : namespace matrixtools { 36 : 37 : class CovarianceMatrix : public ActionShortcut { 38 : public: 39 : static void registerKeywords(Keywords&); 40 : explicit CovarianceMatrix(const ActionOptions&ao); 41 : }; 42 : 43 : PLUMED_REGISTER_ACTION(CovarianceMatrix,"COVARIANCE_MATRIX") 44 : 45 10 : void CovarianceMatrix::registerKeywords(Keywords& keys ) { 46 10 : ActionShortcut::registerKeywords( keys ); 47 20 : keys.add("numbered","ARG","the vectors of data from which we are calculating the covariance"); 48 20 : keys.add("compulsory","WEIGHTS","this keyword takes the label of an action that calculates a vector of values. The elements of this vector " 49 : "are used as weights for the input data points."); 50 20 : keys.addFlag("UNORMALIZED",false,"do not divide by the sum of the weights"); 51 10 : keys.setValueDescription("the covariance matrix"); 52 40 : keys.needsAction("SUM"); keys.needsAction("CUSTOM"); keys.needsAction("VSTACK"); keys.needsAction("TRANSPOSE"); 53 30 : keys.needsAction("ONES"); keys.needsAction("OUTER_PRODUCT"); keys.needsAction("MATRIX_PRODUCT"); 54 10 : } 55 : 56 4 : CovarianceMatrix::CovarianceMatrix(const ActionOptions&ao): 57 : Action(ao), 58 4 : ActionShortcut(ao) 59 : { 60 8 : std::vector<std::string> args; parseVector("ARG",args); 61 4 : unsigned nargs=args.size(); std::string argstr="ARG=" + args[0]; 62 12 : for(unsigned i=1; i<args.size(); ++i) argstr += "," + args[i]; 63 : 64 8 : bool unorm; parseFlag("UNORMALIZED",unorm); std::string wstr; parse("WEIGHTS",wstr); 65 4 : if( !unorm ) { 66 : // Normalize the weights 67 8 : readInputLine( getShortcutLabel() + "_wsum: SUM ARG=" + wstr + " PERIODIC=NO"); 68 8 : readInputLine( getShortcutLabel() + "_weights: CUSTOM ARG=" + wstr + "," + getShortcutLabel() + "_wsum FUNC=x/y PERIODIC=NO"); 69 8 : wstr = getShortcutLabel() + "_weights"; 70 : } 71 : // Make a stack of all the data 72 8 : readInputLine( getShortcutLabel() + "_stack: VSTACK " + argstr ); 73 : // And calculate the covariance matrix by first transposing the stack 74 8 : readInputLine( getShortcutLabel() + "_stackT: TRANSPOSE ARG=" + getShortcutLabel() + "_stack"); 75 : // Create a matrix that holds all the weights 76 4 : std::string str_nargs; Tools::convert( nargs, str_nargs ); 77 8 : readInputLine( getShortcutLabel() + "_ones: ONES SIZE=" + str_nargs ); 78 : // Now create a matrix that holds all the weights 79 8 : readInputLine( getShortcutLabel() + "_matweights: OUTER_PRODUCT ARG=" + getShortcutLabel() + "_ones," + wstr ); 80 : // And multiply the weights by the transpose to get the weighted transpose 81 8 : readInputLine( getShortcutLabel() + "_wT: CUSTOM ARG=" + getShortcutLabel() + "_matweights," + getShortcutLabel() + "_stackT FUNC=x*y PERIODIC=NO"); 82 : // And now calculate the covariance by doing a suitable matrix product 83 8 : readInputLine( getShortcutLabel() + ": MATRIX_PRODUCT ARG=" + getShortcutLabel() + "_wT," + getShortcutLabel() + "_stack"); 84 4 : } 85 : 86 : } 87 : }