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 "core/ActionWithValue.h"
23 : #include "core/ActionWithArguments.h"
24 : #include "core/ActionRegister.h"
25 :
26 : //+PLUMEDOC MCOLVAR CONCATENATE
27 : /*
28 : Join vectors or matrices together
29 :
30 : This action can be used to join sets of scalars, vector or matrices together into a single output value. The following
31 : example shows how this works with vectors and scalars:
32 :
33 : ```plumed
34 : d1: DISTANCE ATOMS=1,2
35 : d2: DISTANCE ATOMS1=3,4 ATOMS2=5,6
36 : v: CONCATENATE ARG=d1,d2
37 : ```
38 :
39 : The output vector here, `v`, has three components. The first of these is the distance between atoms 1 and 2 and the second
40 : and third are the distances between atoms 3 and 4 and 5 and 6 respectively.
41 :
42 : Concatenating two or more vectors is similarly straightfoward. Concatenating matrices is a little more difficult as the
43 : user must explain the axes that each matrix should be joined to. The following input should give some idea as to how
44 : the joining of matrices is done in practise. In this example, we define two groups of atoms and calculate the contact
45 : matrix between all the atoms in the two groups by first calculating three contact matrices that describe the inter and intra
46 : group contacts.
47 :
48 : ```plumed
49 : # Define two groups of atoms
50 : g: GROUP ATOMS=1-5
51 : h: GROUP ATOMS=6-20
52 :
53 : # Calculate the CONTACT_MATRIX for the atoms in group g
54 : cg: CONTACT_MATRIX GROUP=g SWITCH={RATIONAL R_0=0.1}
55 :
56 : # Calculate the CONTACT_MATRIX for the atoms in group h
57 : ch: CONTACT_MATRIX GROUP=h SWITCH={RATIONAL R_0=0.2}
58 :
59 : # Calculate the CONTACT_MATRIX between the atoms in group g and group h
60 : cgh: CONTACT_MATRIX GROUPA=g GROUPB=h SWITCH={RATIONAL R_0=0.15}
61 :
62 : # Now calculate the contact matrix between the atoms in group h and group h
63 : # Notice this is just the transpose of cgh
64 : cghT: TRANSPOSE ARG=cgh
65 :
66 : # And concatenate the matrices together to construct the adjacency matrix between the
67 : # adjacency matrices
68 : m: CONCATENATE ...
69 : MATRIX11=cg MATRIX12=cgh
70 : MATRIX21=cghT MATRIX22=ch
71 : ...
72 : ```
73 :
74 : */
75 : //+ENDPLUMEDOC
76 :
77 : namespace PLMD {
78 : namespace valtools {
79 :
80 : class Concatenate :
81 : public ActionWithValue,
82 : public ActionWithArguments {
83 : private:
84 : bool vectors;
85 : std::vector<unsigned> row_starts;
86 : std::vector<unsigned> col_starts;
87 : public:
88 : static void registerKeywords( Keywords& keys );
89 : /// Constructor
90 : explicit Concatenate(const ActionOptions&);
91 : /// Get the number of derivatives
92 232 : unsigned getNumberOfDerivatives() override {
93 232 : return 0;
94 : }
95 : /// Do the calculation
96 : void calculate() override;
97 : ///
98 : void apply() override;
99 : };
100 :
101 : PLUMED_REGISTER_ACTION(Concatenate,"CONCATENATE")
102 :
103 141 : void Concatenate::registerKeywords( Keywords& keys ) {
104 141 : Action::registerKeywords( keys );
105 141 : ActionWithValue::registerKeywords( keys );
106 141 : ActionWithArguments::registerKeywords( keys );
107 141 : keys.remove("NUMERICAL_DERIVATIVES");
108 282 : keys.addInputKeyword("optional","ARG","scalar/vector","the values that should be concatenated together to form the output vector");
109 282 : keys.addInputKeyword("numbered","MATRIX","scalar/matrix","specify the matrices that you wish to join together into a single matrix");
110 282 : keys.reset_style("MATRIX","compulsory");
111 282 : keys.setValueDescription("vector/matrix","the concatenated vector/matrix that was constructed from the input values");
112 141 : }
113 :
114 69 : Concatenate::Concatenate(const ActionOptions& ao):
115 : Action(ao),
116 : ActionWithValue(ao),
117 69 : ActionWithArguments(ao) {
118 69 : if( getNumberOfArguments()>0 ) {
119 64 : vectors=true;
120 64 : std::vector<std::size_t> shape(1);
121 64 : shape[0]=0;
122 282 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
123 218 : if( getPntrToArgument(i)->getRank()>1 ) {
124 0 : error("cannot concatenate matrix with vectors");
125 : }
126 218 : shape[0] += getPntrToArgument(i)->getNumberOfValues();
127 : }
128 64 : log.printf(" creating vector with %zu elements \n", shape[0] );
129 64 : addValue( shape );
130 64 : bool period=getPntrToArgument(0)->isPeriodic();
131 : std::string min, max;
132 64 : if( period ) {
133 0 : getPntrToArgument(0)->getDomain( min, max );
134 : }
135 218 : for(unsigned i=1; i<getNumberOfArguments(); ++i) {
136 154 : if( period!=getPntrToArgument(i)->isPeriodic() ) {
137 0 : error("periods of input arguments should match");
138 : }
139 154 : if( period ) {
140 : std::string min0, max0;
141 0 : getPntrToArgument(i)->getDomain( min0, max0 );
142 0 : if( min0!=min || max0!=max ) {
143 0 : error("domains of input arguments should match");
144 : }
145 : }
146 : }
147 64 : if( period ) {
148 0 : setPeriodic( min, max );
149 : } else {
150 64 : setNotPeriodic();
151 : }
152 64 : if( getPntrToComponent(0)->getRank()==2 ) {
153 0 : getPntrToComponent(0)->reshapeMatrixStore( shape[1] );
154 : }
155 64 : } else {
156 : unsigned nrows=0, ncols=0;
157 5 : std::vector<Value*> arglist;
158 5 : vectors=false;
159 9 : for(unsigned i=1;; i++) {
160 : unsigned nt_cols=0;
161 : unsigned size_b4 = arglist.size();
162 18 : for(unsigned j=1;; j++) {
163 32 : if( j==10 ) {
164 0 : error("cannot combine more than 9 matrices");
165 : }
166 32 : std::vector<Value*> argn;
167 64 : parseArgumentList("MATRIX", i*10+j, argn);
168 32 : if( argn.size()==0 ) {
169 : break;
170 : }
171 18 : if( argn.size()>1 ) {
172 0 : error("should only be one argument to each matrix keyword");
173 : }
174 18 : if( argn[0]->getRank()!=0 && argn[0]->getRank()!=2 ) {
175 0 : error("input arguments for this action should be matrices");
176 : }
177 18 : arglist.push_back( argn[0] );
178 18 : nt_cols++;
179 18 : if( argn[0]->getRank()==0 ) {
180 0 : log.printf(" %d %d component of composed matrix is scalar labelled %s\n", i, j, argn[0]->getName().c_str() );
181 : } else {
182 18 : log.printf(" %d %d component of composed matrix is %zu by %zu matrix labelled %s\n", i, j, argn[0]->getShape()[0], argn[0]->getShape()[1], argn[0]->getName().c_str() );
183 : }
184 18 : }
185 14 : if( arglist.size()==size_b4 ) {
186 : break;
187 : }
188 9 : if( i==1 ) {
189 : ncols=nt_cols;
190 4 : } else if( nt_cols!=ncols ) {
191 0 : error("should be joining same number of matrices in each row");
192 : }
193 9 : nrows++;
194 9 : }
195 :
196 5 : std::vector<std::size_t> shape(2);
197 5 : shape[0]=0;
198 : unsigned k=0;
199 5 : row_starts.resize( arglist.size() );
200 5 : col_starts.resize( arglist.size() );
201 14 : for(unsigned i=0; i<nrows; ++i) {
202 : unsigned cstart = 0, nr = 1;
203 9 : if( arglist[k]->getRank()==2 ) {
204 9 : nr=arglist[k]->getShape()[0];
205 : }
206 27 : for(unsigned j=0; j<ncols; ++j) {
207 18 : if( arglist[k]->getRank()==0 ) {
208 0 : if( nr!=1 ) {
209 0 : error("mismatched matrix sizes");
210 : }
211 18 : } else if( nrows>1 && arglist[k]->getShape()[0]!=nr ) {
212 0 : error("mismatched matrix sizes");
213 : }
214 18 : row_starts[k] = shape[0];
215 18 : col_starts[k] = cstart;
216 18 : if( arglist[k]->getRank()==0 ) {
217 0 : cstart += 1;
218 : } else {
219 18 : cstart += arglist[k]->getShape()[1];
220 : }
221 18 : k++;
222 : }
223 9 : if( i==0 ) {
224 5 : shape[1]=cstart;
225 4 : } else if( cstart!=shape[1] ) {
226 0 : error("mismatched matrix sizes");
227 : }
228 9 : if( arglist[k-1]->getRank()==0 ) {
229 0 : shape[0] += 1;
230 : } else {
231 9 : shape[0] += arglist[k-1]->getShape()[0];
232 : }
233 : }
234 : // Now request the arguments to make sure we store things we need
235 5 : requestArguments(arglist);
236 5 : addValue( shape );
237 5 : setNotPeriodic();
238 5 : if( getPntrToComponent(0)->getRank()==2 ) {
239 5 : getPntrToComponent(0)->reshapeMatrixStore( shape[1] );
240 : }
241 5 : }
242 69 : }
243 :
244 14823 : void Concatenate::calculate() {
245 14823 : Value* myval = getPntrToComponent(0);
246 14823 : if( vectors ) {
247 : unsigned k=0;
248 71680 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
249 : Value* myarg=getPntrToArgument(i);
250 56926 : unsigned nvals=myarg->getNumberOfValues();
251 518185 : for(unsigned j=0; j<nvals; ++j) {
252 461259 : myval->set( k, myarg->get(j) );
253 461259 : k++;
254 : }
255 : }
256 : } else {
257 : // Retrieve the matrix from input
258 69 : unsigned ncols = myval->getShape()[1];
259 343 : for(unsigned k=0; k<getNumberOfArguments(); ++k) {
260 : Value* argn = getPntrToArgument(k);
261 274 : if( argn->getRank()==0 ) {
262 0 : myval->set( ncols*row_starts[k]+col_starts[k], argn->get() );
263 : } else {
264 274 : std::vector<double> vals;
265 274 : std::vector<std::pair<unsigned,unsigned> > pairs;
266 : bool symmetric=getPntrToArgument(k)->isSymmetric();
267 274 : unsigned nedge=0;
268 274 : getPntrToArgument(k)->retrieveEdgeList( nedge, pairs, vals );
269 11394 : for(unsigned l=0; l<nedge; ++l ) {
270 11120 : unsigned i=pairs[l].first, j=pairs[l].second;
271 11120 : myval->set( ncols*(row_starts[k]+i)+col_starts[k]+j, vals[l] );
272 11120 : if( symmetric ) {
273 2856 : myval->set( ncols*(row_starts[k]+j)+col_starts[k]+i, vals[l] );
274 : }
275 : }
276 274 : }
277 : }
278 : }
279 14823 : }
280 :
281 14816 : void Concatenate::apply() {
282 14816 : if( doNotCalculateDerivatives() || !getPntrToComponent(0)->forcesWereAdded() ) {
283 9216 : return;
284 : }
285 :
286 5600 : Value* val=getPntrToComponent(0);
287 5600 : if( vectors ) {
288 : unsigned k=0;
289 22123 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
290 : Value* myarg=getPntrToArgument(i);
291 16592 : unsigned nvals=myarg->getNumberOfValues();
292 229026 : for(unsigned j=0; j<nvals; ++j) {
293 212434 : myarg->addForce( j, val->getForce(k) );
294 212434 : k++;
295 : }
296 : }
297 : } else {
298 69 : unsigned ncols=val->getShape()[1];
299 343 : for(unsigned k=0; k<getNumberOfArguments(); ++k) {
300 : Value* argn=getPntrToArgument(k);
301 274 : if( argn->getRank()==0 ) {
302 0 : argn->addForce( 0, val->getForce(ncols*row_starts[k]+col_starts[k]) );
303 : } else {
304 274 : unsigned val_ncols=val->getNumberOfColumns();
305 274 : unsigned arg_ncols=argn->getNumberOfColumns();
306 2230 : for(unsigned i=0; i<argn->getShape()[0]; ++i) {
307 : unsigned ncol = argn->getRowLength(i);
308 15996 : for(unsigned j=0; j<ncol; ++j) {
309 14040 : argn->addForce( i*arg_ncols+j, val->getForce( val_ncols*(row_starts[k]+i)+col_starts[k]+argn->getRowIndex(i,j) ), false );
310 : }
311 : }
312 : }
313 : }
314 : }
315 : }
316 :
317 : }
318 : }
|
