Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2020 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/ActionWithVector.h"
23 : #include "core/ParallelTaskManager.h"
24 : #include "core/ActionRegister.h"
25 : #include "tools/ConjugateGradient.h"
26 : #include "tools/SwitchingFunction.h"
27 : #include "tools/OpenMP.h"
28 : #include "tools/Random.h"
29 :
30 : namespace PLMD {
31 : namespace dimred {
32 :
33 : //+PLUMEDOC DIMRED PROJECT_POINTS
34 : /*
35 : Find the projection of a point in a low dimensional space by matching the (transformed) distance between it and a series of reference configurations that were input
36 :
37 : This action and [ARRANGE_POINTS](ARRANGE_POINTS.md) are the workhorses for the implementation of [SKETCHMAP](SKETCHMAP.md) that is provided within PLUMED.
38 : PROJECT_POINTS allows you to provide the low dimensional coordinate $y_\textrm{min}$ at which the following stress function is minimised:
39 :
40 : $$
41 : \chi(y) = \sum_{i=1}^N w_i [ D(X_i, Y) - d(x_i,y) ]^2
42 : $$
43 :
44 : where $Y$ is a set of coordinates in some high dimensional space, $X_i$ is a set of coordinates for one of $N$ landmark points in this high-dimensional space,
45 : $x_i$ is a projection for $X_i$ that we have found in some lower dimensional space and $w_i$ is a weight. The $D$ indicates that we are calculating the dissimilarity between the point $X_i$ and
46 : $Y$, while $d$ represents the distance between the point $x_i$ and $y$. In minimising the expression above we are thus finding the point $y$ at which the distances between $y$ and
47 : each of the projections, $x_i$, of the $N$ landmark points most closely resembles the dissimiarities between $Y$ and the $N$ landmark points in the high-dimensional points.
48 :
49 : The example input below illustrates how you can use PROJECT_POINTS to find the projection of a high-dimensional point in practice.
50 :
51 : ```plumed
52 : # The coordinates of the landmarks in the high dimensional space
53 : d1_ref: CONSTANT VALUES=1.0,2.0,1.5,2.1
54 : d2_ref: CONSTANT VALUES=0.5,0.7,0.2,1.3
55 : d3_ref: CONSTANT VALUES=3.1,2.0,1.5,0.5
56 :
57 : # The weights of the landmark
58 : weights: CONSTANT VALUES=1,1,1,1
59 :
60 : # The projections of the landmarks in the low dimensional space
61 : proj1_ref: CONSTANT VALUES=0.5,0.8,0.2,0.4
62 : proj2_ref: CONSTANT VALUES=0.2,0.9,0.3,0.7
63 :
64 : # Calculate the instantaneous values of the three distances
65 : d1: DISTANCE ATOMS=1,2
66 : d2: DISTANCE ATOMS=3,4
67 : d3: DISTANCE ATOMS=5,6
68 :
69 : # Calculate the distances between the instananeous points and the current positions
70 : ed: EUCLIDEAN_DISTANCE SQUARED ARG2=d1,d2,d3 ARG1=d1_ref,d2_ref,d3_ref
71 :
72 : # And generate the projection
73 : proj: PROJECT_POINTS ARG=proj1_ref,proj2_ref TARGET1=ed WEIGHTS1=weights
74 :
75 : # And output the projection to a file
76 : PRINT ARG=proj.* FILE=colvar
77 : ```
78 :
79 : In this example, we use three distances to define the high dimensional coordinates and have four landmarks points.
80 :
81 : ## Projecting multiple coordinates at once
82 :
83 : The input to the [EUCLIDEAN_DISTANCE](EUCLIDEAN_DISTANCE.md) shortcut in the example in the previous section consisted of three
84 : scalar-valued quantities. The two components output by project points are thus also scalars. By contrast, in the input below four
85 : three-dimenional vectors are input to the [EUCLIDEAN_DISTANCE](EUCLIDEAN_DISTANCE.md) shortcut. The components output by proj and thus
86 : four-dimensional vectors.
87 :
88 : ```plumed
89 : # The coordinates of the landmarks in the high dimensional space
90 : d1_ref: CONSTANT VALUES=1.0,2.0,1.5,2.1
91 : d2_ref: CONSTANT VALUES=0.5,0.7,0.2,1.3
92 : d3_ref: CONSTANT VALUES=3.1,2.0,1.5,0.5
93 :
94 : # The weights of the landmark
95 : weights: CONSTANT VALUES=1,1,1,1
96 :
97 : # The projections of the landmarks in the low dimensional space
98 : proj1_ref: CONSTANT VALUES=0.5,0.8,0.2,0.4
99 : proj2_ref: CONSTANT VALUES=0.2,0.9,0.3,0.7
100 :
101 : # Calculate the instantaneous values of the distances
102 : d1: DISTANCE ATOMS1=1,2 ATOMS2=7,8 ATOMS3=13,14 ATOMS4=19,20
103 : d2: DISTANCE ATOMS1=3,4 ATOMS2=9,10 ATOMS3=15,16 ATOMS4=21,22
104 : d3: DISTANCE ATOMS1=5,6 ATOMS2=11,12 ATOMS3=17,18 ATOMS4=23,24
105 :
106 : # Calculate the distances between the instananeous points and the current positions
107 : ed: EUCLIDEAN_DISTANCE SQUARED ARG1=d1,d2,d3 ARG2=d1_ref,d2_ref,d3_ref
108 :
109 : # And generate the projection
110 : proj: PROJECT_POINTS ARG=proj1_ref,proj2_ref TARGET1=ed WEIGHTS1=weights
111 :
112 : # And output the projection to a file
113 : PRINT ARG=proj.* FILE=colvar
114 : ```
115 :
116 : ## Using RMSD distances
117 :
118 : One can use [RMSD](RMSD.md) distances as the dissimilarities rather than distances in some space of arguments as is illustrated below:
119 :
120 : ```plumed
121 : #SETTINGS INPUTFILES=regtest/trajectories/path_msd/allv.pdb
122 :
123 : # This action reads in the landmarks in the high dimensional space and calculates the
124 : # distances from the instantaneous configuration
125 : rmsd: RMSD SQUARED TYPE=OPTIMAL REFERENCE=regtest/trajectories/path_msd/allv.pdb
126 :
127 : # The weights of the landmarks
128 : weights: ONES SIZE=42
129 :
130 : # The projections of the landmarks in the low dimensional space
131 : X: PDB2CONSTANT ARG=X NOARGS REFERENCE=regtest/trajectories/path_msd/allv.pdb
132 : Y: PDB2CONSTANT ARG=Y NOARGS REFERENCE=regtest/trajectories/path_msd/allv.pdb
133 :
134 : # Generate the projection of the instantaneous coordinates
135 : proj: PROJECT_POINTS ARG=X,Y TARGET1=rmsd WEIGHTS1=weights
136 :
137 : # And output the projection to a file
138 : PRINT ARG=proj.* FILE=colvar
139 : ```
140 :
141 : For ths input there are 42 landmark points and dissimilarities are computed by computing the RMSD distance between the 13 atoms in
142 : each of landmark coordinates and the instantaneous positions of those 13 atoms.
143 :
144 : ## Using transformed distances
145 :
146 : In [SKETCHMAP](SKETCHMAP.md) the stress function that is minimised is not the one given above. Instead of seeking to generate a projection,
147 : $y_\textrm{min}$, which is at a point where the distances between it and each projection the landmarks is the same as the dissimilarities between
148 : the high-dimensional coordinate of the point and the high-dimensional landmarks, the dissimilarities and distances are transformed by functions as illustrated below:
149 :
150 : $$
151 : \chi(y) = \sum_{i=1}^N w_i [ F[D(X_i, Y)] - f[d(x_i,y)] ]^2
152 : $$
153 :
154 : The two functions $F$ and $f$ in this expression are usually different as you can see in the input below:
155 :
156 : ```plumed
157 : # The coordinates of the landmarks in the high dimensional space
158 : d1_ref: CONSTANT VALUES=1.0,2.0,1.5,2.1
159 : d2_ref: CONSTANT VALUES=0.5,0.7,0.2,1.3
160 : d3_ref: CONSTANT VALUES=3.1,2.0,1.5,0.5
161 :
162 : # The weights of the landmark
163 : weights: CONSTANT VALUES=1,1,1,1
164 :
165 : # The projections of the landmarks in the low dimensional space
166 : proj1_ref: CONSTANT VALUES=0.5,0.8,0.2,0.4
167 : proj2_ref: CONSTANT VALUES=0.2,0.9,0.3,0.7
168 :
169 : # Calculate the instantaneous values of the distances
170 : d1: DISTANCE ATOMS1=1,2 ATOMS2=7,8 ATOMS3=13,14 ATOMS4=19,20
171 : d2: DISTANCE ATOMS1=3,4 ATOMS2=9,10 ATOMS3=15,16 ATOMS4=21,22
172 : d3: DISTANCE ATOMS1=5,6 ATOMS2=11,12 ATOMS3=17,18 ATOMS4=23,24
173 :
174 : # Calculate the distances between the instananeous points and the current positions
175 : ed: EUCLIDEAN_DISTANCE SQUARED ARG1=d1,d2,d3 ARG2=d1_ref,d2_ref,d3_ref
176 :
177 : # Transform the dissimilarities by applying the function F
178 : fed: MORE_THAN ARG=ed SQUARED SWITCH={SMAP R_0=4 A=3 B=2}
179 :
180 : # And generate the projection
181 : proj: PROJECT_POINTS ARG=proj1_ref,proj2_ref TARGET1=fed FUNC1={SMAP R_0=4 A=1 B=2} WEIGHTS1=weights
182 :
183 : # And output the projection to a file
184 : PRINT ARG=proj.* FILE=colvar
185 : ```
186 :
187 : In the input above the function, $F$, that is applied on the dissimilarities is implemented using a [MORE_THAN](MORE_THAN.md) action. The function, $f$,
188 : that is applied on the distances in the low-dimensional space is specified using the `FUNC` keyword that is input to PROJECT_POINTS.
189 :
190 : ## Using multiple targets
191 :
192 : At its most complex this action allows you to minimise a stress function such as the one below:
193 :
194 : $$
195 : \chi(y) = \sum_{i=1}^N \sum_{j=1}^M w_{ij} [ F_j[D(X_i, Y)] - f_j[d(x_i,y)] ]^2
196 : $$
197 :
198 : The input below shows how this can be implemented within PLUMED:
199 :
200 : ```plumed
201 : # The coordinates of the landmarks in the high dimensional space
202 : d1_ref: CONSTANT VALUES=1.0,2.0,1.5,2.1
203 : d2_ref: CONSTANT VALUES=0.5,0.7,0.2,1.3
204 : d3_ref: CONSTANT VALUES=3.1,2.0,1.5,0.5
205 :
206 : # The weights of the landmark
207 : weights: CONSTANT VALUES=1,1,1,1
208 : w1: CUSTOM ARG=weights FUNC=0.3*x PERIODIC=NO
209 : w2: CUSTOM ARG=weights FUNC=(1-0.3)*x PERIODIC=NO
210 :
211 : # The projections of the landmarks in the low dimensional space
212 : proj1_ref: CONSTANT VALUES=0.5,0.8,0.2,0.4
213 : proj2_ref: CONSTANT VALUES=0.2,0.9,0.3,0.7
214 :
215 : # Calculate the instantaneous values of the distances
216 : d1: DISTANCE ATOMS1=1,2 ATOMS2=7,8 ATOMS3=13,14 ATOMS4=19,20
217 : d2: DISTANCE ATOMS1=3,4 ATOMS2=9,10 ATOMS3=15,16 ATOMS4=21,22
218 : d3: DISTANCE ATOMS1=5,6 ATOMS2=11,12 ATOMS3=17,18 ATOMS4=23,24
219 :
220 : # Calculate the distances between the instananeous points and the current positions
221 : ed: EUCLIDEAN_DISTANCE SQUARED ARG1=d1,d2,d3 ARG2=d1_ref,d2_ref,d3_ref
222 : # Transform the dissimilarities by applying the function F
223 : fed: MORE_THAN ARG=ed SQUARED SWITCH={SMAP R_0=4 A=3 B=2}
224 :
225 : # And generate the projection
226 : proj: PROJECT_POINTS ...
227 : ARG=proj1_ref,proj2_ref
228 : TARGET1=ed WEIGHTS1=w1 FUNC1={CUSTOM FUNC=1-sqrt(x2) R_0=1.0}
229 : TARGET2=fed WEIGHTS2=w2 FUNC2={SMAP R_0=4 A=1 B=2}
230 : ...
231 :
232 : # And output the projection to a file
233 : PRINT ARG=proj.* FILE=colvar
234 : ```
235 :
236 : Here the sum over $M$ in the expression above has two terms. In the first of these terms $F_1$ is the identity so the
237 : input for `TARGET1` is the output from [EUCLIDEAN_DISTANCE](EUCLIDEAN_DISTANCE.md). $f_1$ is similarly the identity. To
238 : implement the identity here we use the input to `FUNC1` shown above. The input to this function is the input for one of
239 : the switching functions described in the documentation for [LESS_THAN](LESS_THAN.md). What we compute for the transformed
240 : distance is $1-s(d)$ where $s(d)$ is the switching function that is specified in input. Consequently, applying the
241 : function `1-sqrt(x2)` returns the distance.
242 :
243 : The second term in our sum over $M$ in the input above has the dissimilarities and distances transformed by the functions that
244 : we introduced in the previous section.
245 :
246 : */
247 : //+ENDPLUMEDOC
248 :
249 : class ProjectPoints;
250 :
251 : class ProjectPointsInput {
252 : public:
253 : double cgtol;
254 : ProjectPoints* action;
255 : };
256 :
257 : class ProjectPoints : public ActionWithVector {
258 : public:
259 : using input_type = ProjectPointsInput;
260 : using PTM = ParallelTaskManager<ProjectPoints>;
261 : typedef typename PTM::ParallelActionsInput ParallelActionsInput;
262 : typedef typename PTM::ParallelActionsOutput ParallelActionsOutput;
263 : private:
264 : unsigned dimout;
265 : mutable std::vector<unsigned> rowstart;
266 : std::vector<SwitchingFunction> switchingFunction;
267 : ConjugateGradient<ProjectPoints> myminimiser;
268 : PTM taskmanager;
269 : public:
270 : static void registerKeywords( Keywords& keys );
271 : ProjectPoints( const ActionOptions& );
272 0 : unsigned getNumberOfDerivatives() override {
273 0 : return 0;
274 : }
275 : void prepare() override ;
276 : static void performTask( std::size_t task_index, const ProjectPointsInput& actiondata, ParallelActionsInput& input, ParallelActionsOutput& output );
277 : double calculateStress( const std::vector<double>& pp, std::vector<double>& der );
278 : void calculate() override ;
279 168 : void apply() override {}
280 : };
281 :
282 : PLUMED_REGISTER_ACTION(ProjectPoints,"PROJECT_POINTS")
283 :
284 8 : void ProjectPoints::registerKeywords( Keywords& keys ) {
285 8 : ActionWithVector::registerKeywords( keys );
286 16 : keys.addInputKeyword("compulsory","ARG","vector","the projections of the landmark points");
287 16 : keys.addInputKeyword("numbered","TARGET","vector/matrix","the matrix of target quantities that you would like to match");
288 8 : keys.add("numbered","FUNC","a function that is applied on the distances between the points in the low dimensional space");
289 16 : keys.addInputKeyword("numbered","WEIGHTS","vector","the matrix with the weights of the target quantities");
290 8 : keys.add("compulsory","CGTOL","1E-6","the tolerance for the conjugate gradient minimization");
291 16 : keys.addOutputComponent("coord","default","scalar/vector","the coordinates of the points in the low dimensional space");
292 8 : PTM::registerKeywords( keys );
293 8 : }
294 :
295 :
296 2 : ProjectPoints::ProjectPoints( const ActionOptions& ao ) :
297 : Action(ao),
298 : ActionWithVector(ao),
299 2 : rowstart(OpenMP::getNumThreads()),
300 : myminimiser(this),
301 4 : taskmanager(this) {
302 2 : dimout = getNumberOfArguments();
303 2 : unsigned nvals=getPntrToArgument(0)->getNumberOfValues();
304 6 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
305 4 : if( nvals!=getPntrToArgument(i)->getNumberOfValues() ) {
306 0 : error("mismatch between numbers of projections");
307 : }
308 : }
309 2 : std::vector<Value*> args( getArguments() ), target, weights;
310 : std::string sfd, errors;
311 : unsigned ntoproj=0;
312 : // Read in target "distances" and target weights
313 2 : for(unsigned i=1;; ++i) {
314 4 : target.resize(0);
315 8 : if( !parseArgumentList("TARGET",i,target) ) {
316 : break;
317 : }
318 : std::string inum;
319 2 : Tools::convert( i, inum );
320 2 : if( target.size()!=1 ) {
321 0 : error("should only be one value in input to TARGET" + inum );
322 : }
323 2 : if( (target[0]->getRank()!=1 && target[0]->getRank()!=2) || target[0]->hasDerivatives() ) {
324 0 : error("input to TARGET" + inum + " keyword should be a vector/matrix");
325 : }
326 2 : if( target[0]->getShape()[0]!=nvals ) {
327 0 : error("number of rows in target matrix should match number of input coordinates");
328 : }
329 2 : if( i==1 && target[0]->getRank()==1 ) {
330 : ntoproj = 1;
331 1 : } else if( ntoproj==1 && target[0]->getRank()!=1 ) {
332 0 : error("mismatch between numbers of target distances");
333 1 : } else if( i==1 ) {
334 1 : ntoproj = target[0]->getShape()[1];
335 0 : } else if( target[0]->getRank()>1 && ntoproj!=target[0]->getShape()[1] ) {
336 0 : error("mismatch between numbers of target distances");
337 : }
338 4 : if( !parseArgumentList("WEIGHTS",i,weights) ) {
339 0 : error("missing WEIGHTS" + inum + " keyword in input");
340 : }
341 2 : if( weights.size()!=1 ) {
342 0 : error("should only be one value in input to WEIGHTS" + inum );
343 : }
344 2 : if( weights[0]->getRank()!=1 || weights[0]->hasDerivatives() ) {
345 0 : error("input to WEIGHTS" + inum + " keyword should be a vector");
346 : }
347 2 : if( weights[0]->getShape()[0]!=nvals ) {
348 0 : error("number of weights should match number of input coordinates");
349 : }
350 2 : args.push_back( target[0] );
351 2 : args.push_back( weights[0] );
352 2 : bool has_sf = parseNumbered("FUNC",i,sfd);
353 4 : switchingFunction.push_back( SwitchingFunction() );
354 2 : if( !has_sf ) {
355 0 : switchingFunction[i-1].set( "CUSTOM FUNC=1-sqrt(x2) R_0=1.0", errors );
356 : } else {
357 2 : switchingFunction[i-1].set( sfd, errors );
358 2 : if( errors.length()!=0 ) {
359 0 : error("problem reading switching function description " + errors);
360 : }
361 : }
362 2 : log.printf(" %sth term seeks to match transformed distances with those in matrix %s \n", inum.c_str(), target[0]->getName().c_str() );
363 2 : log.printf(" in %sth term distances are transformed by 1-switching function with r_0=%s \n", inum.c_str(), switchingFunction[i-1].description().c_str() );
364 2 : log.printf(" in %sth term weights of matrix elements in stress function are given by %s \n", inum.c_str(), weights[0]->getName().c_str() );
365 2 : }
366 2 : std::vector<std::size_t> shape(1);
367 2 : shape[0]=ntoproj;
368 2 : if( ntoproj==1 ) {
369 1 : shape.resize(0);
370 : }
371 6 : for(unsigned i=0; i<dimout; ++i) {
372 : std::string num;
373 4 : Tools::convert( i+1, num );
374 4 : addComponent( "coord-" + num, shape );
375 8 : componentIsNotPeriodic( "coord-" + num );
376 : }
377 : // Create a list of tasks to perform
378 : double cgtol;
379 2 : parse("CGTOL",cgtol);
380 2 : log.printf(" tolerance for conjugate gradient algorithm equals %f \n",cgtol);
381 2 : requestArguments( args );
382 2 : checkRead();
383 :
384 : // Setup parallel task manager
385 : ProjectPointsInput input;
386 2 : input.cgtol=cgtol;
387 2 : input.action=this;
388 2 : if( ntoproj!=1 ) {
389 1 : taskmanager.setupParallelTaskManager( 0, 0 );
390 : }
391 2 : taskmanager.setActionInput( input );
392 4 : }
393 :
394 169 : void ProjectPoints::prepare() {
395 169 : if( getPntrToComponent(0)->getRank()==0 ) {
396 168 : return;
397 : }
398 :
399 1 : std::vector<std::size_t> shape(1);
400 1 : shape[0] = getPntrToArgument(dimout)->getShape()[0];
401 3 : for(unsigned i=0; i<dimout; ++i) {
402 2 : if( getPntrToComponent(i)->getShape()[0]!=shape[0] ) {
403 2 : getPntrToComponent(i)->setShape(shape);
404 : }
405 : }
406 1 : }
407 :
408 24180 : double ProjectPoints::calculateStress( const std::vector<double>& pp, std::vector<double>& der ) {
409 24180 : unsigned nmatrices = ( getNumberOfArguments() - dimout ) / 2;
410 : double stress=0;
411 :
412 24180 : unsigned t=OpenMP::getThreadNum();
413 24180 : std::vector<double> dtmp( pp.size() );
414 24180 : unsigned nland = getPntrToArgument(0)->getShape()[0];
415 4889418 : for(unsigned i=0; i<nland; ++i) {
416 : // Calculate distance in low dimensional space
417 : double dd2 = 0;
418 14595714 : for(unsigned k=0; k<pp.size(); ++k) {
419 9730476 : dtmp[k] = pp[k] - getPntrToArgument(k)->get(i);
420 9730476 : dd2 += dtmp[k]*dtmp[k];
421 : }
422 :
423 9730476 : for(unsigned k=0; k<nmatrices; ++k ) {
424 : // Now do transformations and calculate differences
425 4865238 : double df, fd = 1. - switchingFunction[k].calculateSqr( dd2, df );
426 : // Get the weight for this connection
427 4865238 : double weight = getPntrToArgument( dimout + 2*k + 1 )->get( i );
428 : // Get the difference for the connection
429 4865238 : double fdiff = fd - getPntrToArgument( dimout + 2*k )->get( rowstart[t]+i );
430 : // Calculate derivatives
431 4865238 : double pref = -2.*weight*fdiff*df;
432 14595714 : for(unsigned n=0; n<pp.size(); ++n) {
433 9730476 : der[n]+=pref*dtmp[n];
434 : }
435 : // Accumulate the total stress
436 4865238 : stress += weight*fdiff*fdiff;
437 : }
438 : }
439 24180 : return stress;
440 : }
441 :
442 668 : void ProjectPoints::performTask( std::size_t task_index, const ProjectPointsInput& actiondata, ParallelActionsInput& input, ParallelActionsOutput& output ) {
443 : // I doubt we are ever going to implement this on the GPU so I think we can leave this declaration here
444 668 : std::vector<double> point( input.ncomponents );
445 668 : std::size_t nland = input.shapedata[0];
446 : std::size_t base = task_index;
447 668 : if( input.ranks[input.ncomponents]==2 ) {
448 500 : auto myargh = ArgumentBookeepingHolder::create( input.ncomponents, input );
449 500 : base = task_index*myargh.shape[1];
450 : }
451 : unsigned closest=0;
452 668 : double mindist = input.inputdata[input.argstarts[input.ncomponents] + base];
453 139112 : for(unsigned i=1; i<nland; ++i) {
454 138444 : double dist = input.inputdata[input.argstarts[input.ncomponents] + base+i];
455 138444 : if( dist<mindist ) {
456 : mindist=dist;
457 : closest=i;
458 : }
459 : }
460 : // Put the initial guess near to the closest landmark -- may wish to use grid here again Sandip??
461 668 : Random random;
462 668 : random.setSeed(-1234);
463 2004 : for(unsigned j=0; j<input.ncomponents; ++j) {
464 1336 : point[j] = input.inputdata[input.argstarts[j] + closest] + (random.RandU01() - 0.5)*0.01;
465 : }
466 : // And do the optimisation
467 668 : actiondata.action->rowstart[OpenMP::getThreadNum()]=task_index;
468 668 : if( input.ranks[input.ncomponents]==2 ) {
469 500 : auto myargh=ArgumentBookeepingHolder::create( input.ncomponents, input );
470 500 : actiondata.action->rowstart[OpenMP::getThreadNum()] = task_index*myargh.shape[1];
471 : }
472 :
473 668 : int code = actiondata.action->myminimiser.minimise( actiondata.cgtol, point, &ProjectPoints::calculateStress );
474 668 : if( code>0 ) {
475 0 : plumed_merror("failure in conjugate gradient minimisation");
476 : }
477 2004 : for(unsigned i=0; i<input.ncomponents; ++i) {
478 1336 : output.values[i] = point[i];
479 : }
480 668 : }
481 :
482 169 : void ProjectPoints::calculate() {
483 169 : if( getPntrToComponent(0)->getRank()==0 ) {
484 : auto myinput =ParallelActionsInput::create( getPbc() );
485 168 : myinput.noderiv = true;
486 168 : myinput.ncomponents = getNumberOfComponents();
487 168 : std::vector<double> input_buffer;
488 168 : getArgumentInputData( input_buffer );
489 168 : myinput.dataSize = input_buffer.size();
490 168 : myinput.inputdata = input_buffer.data();
491 168 : ArgumentsBookkeeping abk;
492 168 : copyArgumentBookeeping( abk );
493 168 : myinput.setupArguments( abk );
494 168 : std::vector<double> buffer;
495 168 : std::vector<double> derivatives, point( getNumberOfComponents() );
496 168 : auto output = ParallelActionsOutput::create( myinput.ncomponents, point.data(), 0, derivatives.data(), 0, buffer.data() );
497 168 : performTask( 0, taskmanager.getActionInput(), myinput, output );
498 504 : for(unsigned i=0; i<point.size(); ++i) {
499 336 : getPntrToComponent(i)->set(point[i]);
500 : }
501 168 : } else {
502 1 : taskmanager.runAllTasks();
503 : }
504 169 : }
505 :
506 : }
507 : }
|
