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 : # Calcuate 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 : # Calcuate 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 instaneous 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 : # Calcuate 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 funciton 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 implemted 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 : # Calcuate 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 funciton 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 : private:
262 : unsigned dimout;
263 : mutable std::vector<unsigned> rowstart;
264 : std::vector<SwitchingFunction> switchingFunction;
265 : ConjugateGradient<ProjectPoints> myminimiser;
266 : PTM taskmanager;
267 : public:
268 : static void registerKeywords( Keywords& keys );
269 : ProjectPoints( const ActionOptions& );
270 0 : unsigned getNumberOfDerivatives() override {
271 0 : return 0;
272 : }
273 : void prepare() override ;
274 : static void performTask( std::size_t task_index, const ProjectPointsInput& actiondata, ParallelActionsInput& input, ParallelActionsOutput& output );
275 : double calculateStress( const std::vector<double>& pp, std::vector<double>& der );
276 : void calculate() override ;
277 168 : void apply() override {}
278 : };
279 :
280 : PLUMED_REGISTER_ACTION(ProjectPoints,"PROJECT_POINTS")
281 :
282 6 : void ProjectPoints::registerKeywords( Keywords& keys ) {
283 6 : ActionWithVector::registerKeywords( keys );
284 12 : keys.addInputKeyword("compulsory","ARG","vector","the projections of the landmark points");
285 12 : keys.addInputKeyword("numbered","TARGET","vector/matrix","the matrix of target quantities that you would like to match");
286 6 : keys.add("numbered","FUNC","a function that is applied on the distances between the points in the low dimensional space");
287 12 : keys.addInputKeyword("numbered","WEIGHTS","vector","the matrix with the weights of the target quantities");
288 6 : keys.add("compulsory","CGTOL","1E-6","the tolerance for the conjugate gradient minimization");
289 12 : keys.addOutputComponent("coord","default","scalar/vector","the coordinates of the points in the low dimensional space");
290 6 : PTM::registerKeywords( keys );
291 6 : }
292 :
293 :
294 2 : ProjectPoints::ProjectPoints( const ActionOptions& ao ) :
295 : Action(ao),
296 : ActionWithVector(ao),
297 2 : rowstart(OpenMP::getNumThreads()),
298 : myminimiser(this),
299 4 : taskmanager(this) {
300 2 : dimout = getNumberOfArguments();
301 2 : unsigned nvals=getPntrToArgument(0)->getNumberOfValues();
302 6 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
303 4 : if( nvals!=getPntrToArgument(i)->getNumberOfValues() ) {
304 0 : error("mismatch between numbers of projections");
305 : }
306 : }
307 2 : std::vector<Value*> args( getArguments() ), target, weights;
308 : std::string sfd, errors;
309 : unsigned ntoproj=0;
310 : // Read in target "distances" and target weights
311 2 : for(unsigned i=1;; ++i) {
312 4 : target.resize(0);
313 8 : if( !parseArgumentList("TARGET",i,target) ) {
314 : break;
315 : }
316 : std::string inum;
317 2 : Tools::convert( i, inum );
318 2 : if( target.size()!=1 ) {
319 0 : error("should only be one value in input to TARGET" + inum );
320 : }
321 2 : if( (target[0]->getRank()!=1 && target[0]->getRank()!=2) || target[0]->hasDerivatives() ) {
322 0 : error("input to TARGET" + inum + " keyword should be a vector/matrix");
323 : }
324 2 : if( target[0]->getShape()[0]!=nvals ) {
325 0 : error("number of rows in target matrix should match number of input coordinates");
326 : }
327 2 : if( i==1 && target[0]->getRank()==1 ) {
328 : ntoproj = 1;
329 1 : } else if( ntoproj==1 && target[0]->getRank()!=1 ) {
330 0 : error("mismatch between numbers of target distances");
331 1 : } else if( i==1 ) {
332 1 : ntoproj = target[0]->getShape()[1];
333 0 : } else if( target[0]->getRank()>1 && ntoproj!=target[0]->getShape()[1] ) {
334 0 : error("mismatch between numbers of target distances");
335 : }
336 4 : if( !parseArgumentList("WEIGHTS",i,weights) ) {
337 0 : error("missing WEIGHTS" + inum + " keyword in input");
338 : }
339 2 : if( weights.size()!=1 ) {
340 0 : error("should only be one value in input to WEIGHTS" + inum );
341 : }
342 2 : if( weights[0]->getRank()!=1 || weights[0]->hasDerivatives() ) {
343 0 : error("input to WEIGHTS" + inum + " keyword should be a vector");
344 : }
345 2 : if( weights[0]->getShape()[0]!=nvals ) {
346 0 : error("number of weights should match number of input coordinates");
347 : }
348 2 : args.push_back( target[0] );
349 2 : args.push_back( weights[0] );
350 2 : bool has_sf = parseNumbered("FUNC",i,sfd);
351 4 : switchingFunction.push_back( SwitchingFunction() );
352 2 : if( !has_sf ) {
353 0 : switchingFunction[i-1].set( "CUSTOM FUNC=1-sqrt(x2) R_0=1.0", errors );
354 : } else {
355 2 : switchingFunction[i-1].set( sfd, errors );
356 2 : if( errors.length()!=0 ) {
357 0 : error("problem reading switching function description " + errors);
358 : }
359 : }
360 2 : log.printf(" %sth term seeks to match tranformed distances with those in matrix %s \n", inum.c_str(), target[0]->getName().c_str() );
361 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() );
362 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() );
363 2 : }
364 2 : std::vector<std::size_t> shape(1);
365 2 : shape[0]=ntoproj;
366 2 : if( ntoproj==1 ) {
367 1 : shape.resize(0);
368 : }
369 6 : for(unsigned i=0; i<dimout; ++i) {
370 : std::string num;
371 4 : Tools::convert( i+1, num );
372 4 : addComponent( "coord-" + num, shape );
373 8 : componentIsNotPeriodic( "coord-" + num );
374 : }
375 : // Create a list of tasks to perform
376 : double cgtol;
377 2 : parse("CGTOL",cgtol);
378 2 : log.printf(" tolerance for conjugate gradient algorithm equals %f \n",cgtol);
379 2 : requestArguments( args );
380 2 : checkRead();
381 :
382 : // Setup parallel task manager
383 : ProjectPointsInput input;
384 2 : input.cgtol=cgtol;
385 2 : input.action=this;
386 2 : if( ntoproj!=1 ) {
387 1 : taskmanager.setupParallelTaskManager( 0, 0 );
388 : }
389 2 : taskmanager.setActionInput( input );
390 2 : }
391 :
392 169 : void ProjectPoints::prepare() {
393 169 : if( getPntrToComponent(0)->getRank()==0 ) {
394 168 : return;
395 : }
396 :
397 1 : std::vector<std::size_t> shape(1);
398 1 : shape[0] = getPntrToArgument(dimout)->getShape()[0];
399 3 : for(unsigned i=0; i<dimout; ++i) {
400 2 : if( getPntrToComponent(i)->getShape()[0]!=shape[0] ) {
401 2 : getPntrToComponent(i)->setShape(shape);
402 : }
403 : }
404 : }
405 :
406 24180 : double ProjectPoints::calculateStress( const std::vector<double>& pp, std::vector<double>& der ) {
407 24180 : unsigned nmatrices = ( getNumberOfArguments() - dimout ) / 2;
408 : double stress=0;
409 :
410 24180 : unsigned t=OpenMP::getThreadNum();
411 24180 : std::vector<double> dtmp( pp.size() );
412 24180 : unsigned nland = getPntrToArgument(0)->getShape()[0];
413 4889418 : for(unsigned i=0; i<nland; ++i) {
414 : // Calculate distance in low dimensional space
415 : double dd2 = 0;
416 14595714 : for(unsigned k=0; k<pp.size(); ++k) {
417 9730476 : dtmp[k] = pp[k] - getPntrToArgument(k)->get(i);
418 9730476 : dd2 += dtmp[k]*dtmp[k];
419 : }
420 :
421 9730476 : for(unsigned k=0; k<nmatrices; ++k ) {
422 : // Now do transformations and calculate differences
423 4865238 : double df, fd = 1. - switchingFunction[k].calculateSqr( dd2, df );
424 : // Get the weight for this connection
425 4865238 : double weight = getPntrToArgument( dimout + 2*k + 1 )->get( i );
426 : // Get the difference for the connection
427 4865238 : double fdiff = fd - getPntrToArgument( dimout + 2*k )->get( rowstart[t]+i );
428 : // Calculate derivatives
429 4865238 : double pref = -2.*weight*fdiff*df;
430 14595714 : for(unsigned n=0; n<pp.size(); ++n) {
431 9730476 : der[n]+=pref*dtmp[n];
432 : }
433 : // Accumulate the total stress
434 4865238 : stress += weight*fdiff*fdiff;
435 : }
436 : }
437 24180 : return stress;
438 : }
439 :
440 668 : void ProjectPoints::performTask( std::size_t task_index, const ProjectPointsInput& actiondata, ParallelActionsInput& input, ParallelActionsOutput& output ) {
441 : // I doubt we are ever going to implement this on the GPU so I think we can leave this declaration here
442 668 : std::vector<double> point( input.ncomponents );
443 668 : std::size_t nland = input.shapedata[0];
444 : std::size_t base = task_index;
445 668 : if( input.ranks[input.ncomponents]==2 ) {
446 500 : auto myargh = ArgumentBookeepingHolder::create( input.ncomponents, input );
447 500 : base = task_index*myargh.shape[1];
448 : }
449 : unsigned closest=0;
450 668 : double mindist = input.inputdata[input.argstarts[input.ncomponents] + base];
451 139112 : for(unsigned i=1; i<nland; ++i) {
452 138444 : double dist = input.inputdata[input.argstarts[input.ncomponents] + base+i];
453 138444 : if( dist<mindist ) {
454 : mindist=dist;
455 : closest=i;
456 : }
457 : }
458 : // Put the initial guess near to the closest landmark -- may wish to use grid here again Sandip??
459 668 : Random random;
460 668 : random.setSeed(-1234);
461 2004 : for(unsigned j=0; j<input.ncomponents; ++j) {
462 1336 : point[j] = input.inputdata[input.argstarts[j] + closest] + (random.RandU01() - 0.5)*0.01;
463 : }
464 : // And do the optimisation
465 668 : actiondata.action->rowstart[OpenMP::getThreadNum()]=task_index;
466 668 : if( input.ranks[input.ncomponents]==2 ) {
467 500 : auto myargh=ArgumentBookeepingHolder::create( input.ncomponents, input );
468 500 : actiondata.action->rowstart[OpenMP::getThreadNum()] = task_index*myargh.shape[1];
469 : }
470 668 : actiondata.action->myminimiser.minimise( actiondata.cgtol, point, &ProjectPoints::calculateStress );
471 2004 : for(unsigned i=0; i<input.ncomponents; ++i) {
472 1336 : output.values[i] = point[i];
473 : }
474 668 : }
475 :
476 169 : void ProjectPoints::calculate() {
477 169 : if( getPntrToComponent(0)->getRank()==0 ) {
478 : auto myinput =ParallelActionsInput::create( getPbc() );
479 168 : myinput.noderiv = true;
480 168 : myinput.ncomponents = getNumberOfComponents();
481 : std::vector<double> input_buffer;
482 168 : getInputData( input_buffer );
483 168 : myinput.dataSize = input_buffer.size();
484 168 : myinput.inputdata = input_buffer.data();
485 168 : ArgumentsBookkeeping abk;
486 168 : abk.setupArguments( this );
487 168 : myinput.setupArguments( abk );
488 : std::vector<double> buffer;
489 168 : std::vector<double> derivatives, point( getNumberOfComponents() );
490 168 : auto output = ParallelActionsOutput::create( myinput.ncomponents, point.data(), 0, derivatives.data(), 0, buffer.data() );
491 168 : performTask( 0, taskmanager.getActionInput(), myinput, output );
492 504 : for(unsigned i=0; i<point.size(); ++i) {
493 336 : getPntrToComponent(i)->set(point[i]);
494 : }
495 168 : } else {
496 1 : taskmanager.runAllTasks();
497 : }
498 169 : }
499 :
500 : }
501 : }
|
