Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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 :
23 : #include "Grid.h"
24 : #include "Tools.h"
25 : #include "core/Value.h"
26 : #include "File.h"
27 : #include "Exception.h"
28 : #include "KernelFunctions.h"
29 : #include "RootFindingBase.h"
30 : #include "Communicator.h"
31 :
32 : #include <vector>
33 : #include <cmath>
34 : #include <iostream>
35 : #include <sstream>
36 : #include <cstdio>
37 : #include <cfloat>
38 :
39 : using namespace std;
40 : namespace PLMD {
41 :
42 1102 : Grid::Grid(const std::string& funcl, const std::vector<Value*> & args, const vector<std::string> & gmin,
43 5510 : const vector<std::string> & gmax, const vector<unsigned> & nbin, bool dospline, bool usederiv, bool doclear) {
44 : // various checks
45 1102 : plumed_massert(args.size()==gmin.size(),"grid min dimensions in input do not match number of arguments");
46 1102 : plumed_massert(args.size()==nbin.size(),"number of bins on input do not match number of arguments");
47 1102 : plumed_massert(args.size()==gmax.size(),"grid max dimensions in input do not match number of arguments");
48 1102 : unsigned dim=gmax.size();
49 1102 : std::vector<std::string> names;
50 : std::vector<bool> isperiodic;
51 1102 : std::vector<string> pmin,pmax;
52 1102 : names.resize( dim );
53 1102 : isperiodic.resize( dim );
54 1102 : pmin.resize( dim );
55 1102 : pmax.resize( dim );
56 3886 : for(unsigned int i=0; i<dim; ++i) {
57 2784 : names[i]=args[i]->getName();
58 1392 : if( args[i]->isPeriodic() ) {
59 : isperiodic[i]=true;
60 914 : args[i]->getDomain( pmin[i], pmax[i] );
61 : } else {
62 : isperiodic[i]=false;
63 : pmin[i]="0.";
64 : pmax[i]="0.";
65 : }
66 : }
67 : // this is a value-independent initializator
68 1102 : Init(funcl,names,gmin,gmax,nbin,dospline,usederiv,doclear,isperiodic,pmin,pmax);
69 1102 : }
70 :
71 66 : Grid::Grid(const std::string& funcl, const std::vector<string> &names, const std::vector<std::string> & gmin,
72 330 : const vector<std::string> & gmax, const std::vector<unsigned> & nbin, bool dospline, bool usederiv, bool doclear, const std::vector<bool> &isperiodic, const std::vector<std::string> &pmin, const std::vector<std::string> &pmax ) {
73 : // this calls the initializator
74 66 : Init(funcl,names,gmin,gmax,nbin,dospline,usederiv,doclear,isperiodic,pmin,pmax);
75 66 : }
76 :
77 1168 : void Grid::Init(const std::string& funcl, const std::vector<std::string> &names, const vector<std::string> & gmin,
78 : const std::vector<std::string> & gmax, const std::vector<unsigned> & nbin, bool dospline, bool usederiv, bool doclear,
79 : const std::vector<bool> &isperiodic, const std::vector<std::string> &pmin, const std::vector<std::string> &pmax ) {
80 1168 : contour_location=0.0; fmt_="%14.9f";
81 : // various checks
82 1168 : plumed_massert(names.size()==gmin.size(),"grid dimensions in input do not match number of arguments");
83 1168 : plumed_massert(names.size()==nbin.size(),"grid dimensions in input do not match number of arguments");
84 1168 : plumed_massert(names.size()==gmax.size(),"grid dimensions in input do not match number of arguments");
85 1168 : dimension_=gmax.size();
86 1168 : str_min_=gmin; str_max_=gmax;
87 1168 : argnames.resize( dimension_ );
88 1168 : min_.resize( dimension_ );
89 1168 : max_.resize( dimension_ );
90 1168 : pbc_.resize( dimension_ );
91 4084 : for(unsigned int i=0; i<dimension_; ++i) {
92 1458 : argnames[i]=names[i];
93 1458 : if( isperiodic[i] ) {
94 : pbc_[i]=true;
95 : str_min_[i]=pmin[i];
96 : str_max_[i]=pmax[i];
97 : } else {
98 : pbc_[i]=false;
99 : }
100 1458 : Tools::convert(str_min_[i],min_[i]);
101 1458 : Tools::convert(str_max_[i],max_[i]);
102 1458 : funcname=funcl;
103 2916 : plumed_massert(max_[i]>min_[i],"maximum in grid must be larger than minimum");
104 1458 : plumed_massert(nbin[i]>0,"number of grid points must be greater than zero");
105 : }
106 1168 : nbin_=nbin;
107 1168 : dospline_=dospline;
108 1168 : usederiv_=usederiv;
109 1168 : if(dospline_) plumed_assert(dospline_==usederiv_);
110 1168 : maxsize_=1;
111 4084 : for(unsigned int i=0; i<dimension_; ++i) {
112 7290 : dx_.push_back( (max_[i]-min_[i])/static_cast<double>( nbin_[i] ) );
113 4269 : if( !pbc_[i] ) { max_[i] += dx_[i]; nbin_[i] += 1; }
114 1458 : maxsize_*=nbin_[i];
115 : }
116 1168 : if(doclear) clear();
117 1168 : }
118 :
119 1157 : void Grid::clear() {
120 1157 : grid_.resize(maxsize_);
121 1157 : if(usederiv_) der_.resize(maxsize_);
122 5396383 : for(index_t i=0; i<maxsize_; ++i) {
123 2697613 : grid_[i]=0.0;
124 2697613 : if(usederiv_) {
125 568867 : (der_[i]).resize(dimension_);
126 2825727 : for(unsigned int j=0; j<dimension_; ++j) der_[i][j]=0.0;
127 : }
128 : }
129 1157 : }
130 :
131 96710 : vector<std::string> Grid::getMin() const {
132 96710 : return str_min_;
133 : }
134 :
135 513 : vector<std::string> Grid::getMax() const {
136 513 : return str_max_;
137 : }
138 :
139 7270717 : vector<double> Grid::getDx() const {
140 7270717 : return dx_;
141 : }
142 :
143 28 : double Grid::getBinVolume() const {
144 : double vol=1.;
145 224 : for(unsigned i=0; i<dx_.size(); ++i) vol*=dx_[i];
146 28 : return vol;
147 : }
148 :
149 1933 : vector<bool> Grid::getIsPeriodic() const {
150 1933 : return pbc_;
151 : }
152 :
153 690334 : vector<unsigned> Grid::getNbin() const {
154 690334 : return nbin_;
155 : }
156 :
157 7194 : vector<string> Grid::getArgNames() const {
158 7194 : return argnames;
159 : }
160 :
161 :
162 16449094 : Grid::index_t Grid::getSize() const {
163 16449094 : return maxsize_;
164 : }
165 :
166 15184774 : unsigned Grid::getDimension() const {
167 15184774 : return dimension_;
168 : }
169 :
170 : // we are flattening arrays using a column-major order
171 5810807 : Grid::index_t Grid::getIndex(const vector<unsigned> & indices) const {
172 : plumed_dbg_assert(indices.size()==dimension_);
173 29028899 : for(unsigned int i=0; i<dimension_; i++)
174 34827138 : if(indices[i]>=nbin_[i]) {
175 : std::string is;
176 0 : Tools::convert(i,is);
177 0 : std::string msg="ERROR: the system is looking for a value outside the grid along the " + is + " ("+getArgNames()[i]+")";
178 0 : plumed_merror(msg+" index!");
179 : }
180 11621614 : index_t index=indices[dimension_-1];
181 17407285 : for(unsigned int i=dimension_-1; i>0; --i) {
182 17394717 : index=index*nbin_[i-1]+indices[i-1];
183 : }
184 5810807 : return index;
185 : }
186 :
187 94405 : Grid::index_t Grid::getIndex(const vector<double> & x) const {
188 : plumed_dbg_assert(x.size()==dimension_);
189 188810 : return getIndex(getIndices(x));
190 : }
191 :
192 : // we are flattening arrays using a column-major order
193 17990119 : vector<unsigned> Grid::getIndices(index_t index) const {
194 17990119 : vector<unsigned> indices(dimension_);
195 : index_t kk=index;
196 35980238 : indices[0]=(index%nbin_[0]);
197 17990119 : for(unsigned int i=1; i<dimension_-1; ++i) {
198 0 : kk=(kk-indices[i-1])/nbin_[i-1];
199 0 : indices[i]=(kk%nbin_[i]);
200 : }
201 17990119 : if(dimension_>=2) {
202 69601132 : indices[dimension_-1]=((kk-indices[dimension_-2])/nbin_[dimension_-2]);
203 : }
204 17990119 : return indices;
205 : }
206 :
207 1904931 : vector<unsigned> Grid::getIndices(const vector<double> & x) const {
208 : plumed_dbg_assert(x.size()==dimension_);
209 : vector<unsigned> indices;
210 9511767 : for(unsigned int i=0; i<dimension_; ++i) {
211 19017090 : indices.push_back(unsigned(floor((x[i]-min_[i])/dx_[i])));
212 : }
213 1904931 : return indices;
214 : }
215 :
216 6813565 : vector<double> Grid::getPoint(const vector<unsigned> & indices) const {
217 : plumed_dbg_assert(indices.size()==dimension_);
218 : vector<double> x;
219 33360957 : for(unsigned int i=0; i<dimension_; ++i) {
220 66368480 : x.push_back(min_[i]+(double)(indices[i])*dx_[i]);
221 : }
222 6813565 : return x;
223 : }
224 :
225 5717855 : vector<double> Grid::getPoint(index_t index) const {
226 : plumed_dbg_assert(index<maxsize_);
227 11435710 : return getPoint(getIndices(index));
228 : }
229 :
230 903316 : vector<double> Grid::getPoint(const vector<double> & x) const {
231 : plumed_dbg_assert(x.size()==dimension_);
232 1806632 : return getPoint(getIndices(x));
233 : }
234 :
235 1101020 : void Grid::getPoint(index_t index,std::vector<double> & point) const {
236 : plumed_dbg_assert(index<maxsize_);
237 2202040 : getPoint(getIndices(index),point);
238 1101020 : }
239 :
240 1101020 : void Grid::getPoint(const std::vector<unsigned> & indices,std::vector<double> & point) const {
241 : plumed_dbg_assert(indices.size()==dimension_);
242 : plumed_dbg_assert(point.size()==dimension_);
243 5499916 : for(unsigned int i=0; i<dimension_; ++i) {
244 10997240 : point[i]=(min_[i]+(double)(indices[i])*dx_[i]);
245 : }
246 1101020 : }
247 :
248 0 : void Grid::getPoint(const std::vector<double> & x,std::vector<double> & point) const {
249 : plumed_dbg_assert(x.size()==dimension_);
250 0 : getPoint(getIndices(x),point);
251 0 : }
252 :
253 23046 : vector<Grid::index_t> Grid::getNeighbors
254 : (const vector<unsigned> &indices,const vector<unsigned> &nneigh)const {
255 : plumed_dbg_assert(indices.size()==dimension_ && nneigh.size()==dimension_);
256 :
257 : vector<index_t> neighbors;
258 23046 : vector<unsigned> small_bin(dimension_);
259 :
260 : unsigned small_nbin=1;
261 114830 : for(unsigned j=0; j<dimension_; ++j) {
262 137676 : small_bin[j]=(2*nneigh[j]+1);
263 45892 : small_nbin*=small_bin[j];
264 : }
265 :
266 23046 : vector<unsigned> small_indices(dimension_);
267 : vector<unsigned> tmp_indices;
268 2492454 : for(unsigned index=0; index<small_nbin; ++index) {
269 1234704 : tmp_indices.resize(dimension_);
270 : unsigned kk=index;
271 2469408 : small_indices[0]=(index%small_bin[0]);
272 1234704 : for(unsigned i=1; i<dimension_-1; ++i) {
273 0 : kk=(kk-small_indices[i-1])/small_bin[i-1];
274 0 : small_indices[i]=(kk%small_bin[i]);
275 : }
276 1234704 : if(dimension_>=2) {
277 4927600 : small_indices[dimension_-1]=((kk-small_indices[dimension_-2])/small_bin[dimension_-2]);
278 : }
279 : unsigned ll=0;
280 6167912 : for(unsigned i=0; i<dimension_; ++i) {
281 9866416 : int i0=small_indices[i]-nneigh[i]+indices[i];
282 2466604 : if(!pbc_[i] && i0<0) continue;
283 2580690 : if(!pbc_[i] && i0>=static_cast<int>(nbin_[i])) continue;
284 4822879 : if( pbc_[i] && i0<0) i0=nbin_[i]-(-i0)%nbin_[i];
285 4816932 : if( pbc_[i] && i0>=static_cast<int>(nbin_[i])) i0%=nbin_[i];
286 4929516 : tmp_indices[ll]=static_cast<unsigned>(i0);
287 2464758 : ll++;
288 : }
289 1234704 : tmp_indices.resize(ll);
290 2467562 : if(tmp_indices.size()==dimension_) {neighbors.push_back(getIndex(tmp_indices));}
291 : }
292 23046 : return neighbors;
293 : }
294 :
295 3894 : vector<Grid::index_t> Grid::getNeighbors
296 : (const vector<double> & x,const vector<unsigned> & nneigh)const {
297 : plumed_dbg_assert(x.size()==dimension_ && nneigh.size()==dimension_);
298 7788 : return getNeighbors(getIndices(x),nneigh);
299 : }
300 :
301 19152 : vector<Grid::index_t> Grid::getNeighbors
302 : (index_t index,const vector<unsigned> & nneigh)const {
303 : plumed_dbg_assert(index<maxsize_ && nneigh.size()==dimension_);
304 38304 : return getNeighbors(getIndices(index),nneigh);
305 : }
306 :
307 903316 : vector<Grid::index_t> Grid::getSplineNeighbors(const vector<unsigned> & indices)const {
308 : plumed_dbg_assert(indices.size()==dimension_);
309 : vector<index_t> neighbors;
310 1806632 : unsigned nneigh=unsigned(pow(2.0,int(dimension_)));
311 :
312 8122584 : for(unsigned int i=0; i<nneigh; ++i) {
313 : unsigned tmp=i;
314 : vector<unsigned> nindices;
315 18040910 : for(unsigned int j=0; j<dimension_; ++j) {
316 14431276 : unsigned i0=tmp%2+indices[j];
317 7215638 : tmp/=2;
318 7226072 : if(!pbc_[j] && i0==nbin_[j]) continue;
319 14420836 : if( pbc_[j] && i0==nbin_[j]) i0=0;
320 7215632 : nindices.push_back(i0);
321 : }
322 7219262 : if(nindices.size()==dimension_) neighbors.push_back(getIndex(nindices));
323 : }
324 903316 : return neighbors;
325 : }
326 :
327 9576 : vector<Grid::index_t> Grid::getNearestNeighbors(const index_t index) const {
328 : vector<index_t> nearest_neighs = vector<index_t>();
329 47880 : for (unsigned i = 0; i < dimension_; i++) {
330 19152 : vector<unsigned> neighsneeded = vector<unsigned>(dimension_, 0);
331 38304 : neighsneeded[i] = 1;
332 19152 : vector<index_t> singledim_nearest_neighs = getNeighbors(index, neighsneeded);
333 205770 : for (unsigned j = 0; j < singledim_nearest_neighs.size(); j++) {
334 55822 : index_t neigh = singledim_nearest_neighs[j];
335 55822 : if (neigh != index) {
336 36670 : nearest_neighs.push_back(neigh);
337 : }
338 : }
339 : }
340 9576 : return nearest_neighs;
341 : }
342 :
343 0 : vector<Grid::index_t> Grid::getNearestNeighbors(const vector<unsigned> &indices) const {
344 : plumed_dbg_assert(indices.size() == dimension_);
345 0 : return getNearestNeighbors(getIndex(indices));
346 : }
347 :
348 :
349 0 : void Grid::addKernel( const KernelFunctions& kernel ) {
350 : plumed_dbg_assert( kernel.ndim()==dimension_ );
351 0 : std::vector<unsigned> nneighb=kernel.getSupport( dx_ );
352 0 : std::vector<index_t> neighbors=getNeighbors( kernel.getCenter(), nneighb );
353 0 : std::vector<double> xx( dimension_ ); std::vector<Value*> vv( dimension_ );
354 : std::string str_min, str_max;
355 0 : for(unsigned i=0; i<dimension_; ++i) {
356 0 : vv[i]=new Value();
357 0 : if( pbc_[i] ) {
358 0 : Tools::convert(min_[i],str_min);
359 0 : Tools::convert(max_[i],str_max);
360 0 : vv[i]->setDomain( str_min, str_max );
361 : } else {
362 0 : vv[i]->setNotPeriodic();
363 : }
364 : }
365 :
366 0 : std::vector<double> der( dimension_ );
367 0 : for(unsigned i=0; i<neighbors.size(); ++i) {
368 0 : index_t ineigh=neighbors[i];
369 0 : getPoint( ineigh, xx );
370 0 : for(unsigned j=0; j<dimension_; ++j) vv[j]->set(xx[j]);
371 0 : double newval = kernel.evaluate( vv, der, usederiv_ );
372 0 : if( usederiv_ ) addValueAndDerivatives( ineigh, newval, der );
373 0 : else addValue( ineigh, newval );
374 : }
375 :
376 0 : for(unsigned i=0; i<dimension_; ++i) delete vv[i];
377 0 : }
378 :
379 11458578 : double Grid::getValue(index_t index) const {
380 : plumed_dbg_assert(index<maxsize_);
381 11458578 : return grid_[index];
382 : }
383 :
384 770 : double Grid::getMinValue() const {
385 : double minval;
386 : minval=DBL_MAX;
387 4365888 : for(index_t i=0; i<grid_.size(); ++i) {
388 2182559 : if(grid_[i]<minval)minval=grid_[i];
389 : }
390 770 : return minval;
391 : }
392 :
393 494 : double Grid::getMaxValue() const {
394 : double maxval;
395 : maxval=DBL_MIN;
396 2924486 : for(index_t i=0; i<grid_.size(); ++i) {
397 1461996 : if(grid_[i]>maxval)maxval=grid_[i];
398 : }
399 494 : return maxval;
400 : }
401 :
402 :
403 873916 : double Grid::getValue(const vector<unsigned> & indices) const {
404 873916 : return getValue(getIndex(indices));
405 : }
406 :
407 2264 : double Grid::getValue(const vector<double> & x) const {
408 2264 : if(!dospline_) {
409 18 : return getValue(getIndex(x));
410 : } else {
411 2246 : vector<double> der(dimension_);
412 2246 : return getValueAndDerivatives(x,der);
413 : }
414 : }
415 :
416 4146302 : double Grid::getValueAndDerivatives
417 : (index_t index, vector<double>& der) const {
418 : plumed_dbg_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
419 4146302 : der=der_[index];
420 4146302 : return grid_[index];
421 : }
422 :
423 0 : double Grid::getValueAndDerivatives
424 : (const vector<unsigned> & indices, vector<double>& der) const {
425 0 : return getValueAndDerivatives(getIndex(indices),der);
426 : }
427 :
428 903316 : double Grid::getValueAndDerivatives
429 : (const vector<double> & x, vector<double>& der) const {
430 : plumed_dbg_assert(der.size()==dimension_ && usederiv_);
431 :
432 903316 : if(dospline_) {
433 : double X,X2,X3,value;
434 903316 : vector<double> fd(dimension_);
435 903316 : vector<double> C(dimension_);
436 903316 : vector<double> D(dimension_);
437 903316 : vector<double> dder(dimension_);
438 : // reset
439 : value=0.0;
440 4512950 : for(unsigned int i=0; i<dimension_; ++i) der[i]=0.0;
441 :
442 903316 : vector<unsigned> indices=getIndices(x);
443 903316 : vector<index_t> neigh=getSplineNeighbors(indices);
444 903316 : vector<double> xfloor=getPoint(x);
445 :
446 : // loop over neighbors
447 12635516 : for(unsigned int ipoint=0; ipoint<neigh.size(); ++ipoint) {
448 7219256 : double grid=getValueAndDerivatives(neigh[ipoint],dder);
449 3609628 : vector<unsigned> nindices=getIndices(neigh[ipoint]);
450 : double ff=1.0;
451 :
452 18040892 : for(unsigned j=0; j<dimension_; ++j) {
453 : int x0=1;
454 21646896 : if(nindices[j]==indices[j]) x0=0;
455 21646896 : double dx=getDx()[j];
456 14431264 : X=fabs((x[j]-xfloor[j])/dx-(double)x0);
457 7215632 : X2=X*X;
458 7215632 : X3=X2*X;
459 : double yy;
460 7215632 : if(fabs(grid)<0.0000001) yy=0.0;
461 1004332 : else yy=-dder[j]/grid;
462 14431264 : C[j]=(1.0-3.0*X2+2.0*X3) - (x0?-1.0:1.0)*yy*(X-2.0*X2+X3)*dx;
463 14431264 : D[j]=( -6.0*X +6.0*X2) - (x0?-1.0:1.0)*yy*(1.0-4.0*X +3.0*X2)*dx;
464 7215632 : D[j]*=(x0?-1.0:1.0)/dx;
465 7215632 : ff*=C[j];
466 : }
467 18040892 : for(unsigned j=0; j<dimension_; ++j) {
468 14431264 : fd[j]=D[j];
469 36067288 : for(unsigned i=0; i<dimension_; ++i) if(i!=j) fd[j]*=C[i];
470 : }
471 3609628 : value+=grid*ff;
472 25256524 : for(unsigned j=0; j<dimension_; ++j) der[j]+=grid*fd[j];
473 : }
474 : return value;
475 : } else {
476 0 : return getValueAndDerivatives(getIndex(x),der);
477 : }
478 : }
479 :
480 6066982 : void Grid::setValue(index_t index, double value) {
481 : plumed_dbg_assert(index<maxsize_ && !usederiv_);
482 6066982 : grid_[index]=value;
483 6066982 : }
484 :
485 0 : void Grid::setValue(const vector<unsigned> & indices, double value) {
486 0 : setValue(getIndex(indices),value);
487 0 : }
488 :
489 1507828 : void Grid::setValueAndDerivatives
490 : (index_t index, double value, vector<double>& der) {
491 : plumed_dbg_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
492 1507828 : grid_[index]=value;
493 1507828 : der_[index]=der;
494 1507828 : }
495 :
496 0 : void Grid::setValueAndDerivatives
497 : (const vector<unsigned> & indices, double value, vector<double>& der) {
498 0 : setValueAndDerivatives(getIndex(indices),value,der);
499 0 : }
500 :
501 0 : void Grid::addValue(index_t index, double value) {
502 : plumed_dbg_assert(index<maxsize_ && !usederiv_);
503 0 : grid_[index]+=value;
504 0 : }
505 :
506 0 : void Grid::addValue(const vector<unsigned> & indices, double value) {
507 0 : addValue(getIndex(indices),value);
508 0 : }
509 :
510 1159337 : void Grid::addValueAndDerivatives
511 : (index_t index, double value, vector<double>& der) {
512 : plumed_dbg_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
513 1159337 : grid_[index]+=value;
514 8101628 : for(unsigned int i=0; i<dimension_; ++i) der_[index][i]+=der[i];
515 1159337 : }
516 :
517 0 : void Grid::addValueAndDerivatives
518 : (const vector<unsigned> & indices, double value, vector<double>& der) {
519 0 : addValueAndDerivatives(getIndex(indices),value,der);
520 0 : }
521 :
522 546 : void Grid::scaleAllValuesAndDerivatives( const double& scalef ) {
523 546 : if(usederiv_) {
524 41930 : for(index_t i=0; i<grid_.size(); ++i) {
525 20962 : grid_[i]*=scalef;
526 104810 : for(unsigned j=0; j<dimension_; ++j) der_[i][j]*=scalef;
527 : }
528 : } else {
529 3024124 : for(index_t i=0; i<grid_.size(); ++i) grid_[i]*=scalef;
530 : }
531 546 : }
532 :
533 0 : void Grid::logAllValuesAndDerivatives( const double& scalef ) {
534 0 : if(usederiv_) {
535 0 : for(index_t i=0; i<grid_.size(); ++i) {
536 0 : grid_[i] = scalef*log(grid_[i]);
537 0 : for(unsigned j=0; j<dimension_; ++j) der_[i][j] = scalef/der_[i][j];
538 : }
539 : } else {
540 0 : for(index_t i=0; i<grid_.size(); ++i) grid_[i] = scalef*log(grid_[i]);
541 : }
542 0 : }
543 :
544 1191 : void Grid::setMinToZero() {
545 1191 : double min=grid_[0];
546 6887715 : for(index_t i=1; i<grid_.size(); ++i) if(grid_[i]<min) min=grid_[i];
547 6890097 : for(index_t i=0; i<grid_.size(); ++i) grid_[i] -= min;
548 1191 : }
549 :
550 1 : void Grid::applyFunctionAllValuesAndDerivatives( double (*func)(double val), double (*funcder)(double valder) ) {
551 1 : if(usederiv_) {
552 1801 : for(index_t i=0; i<grid_.size(); ++i) {
553 900 : grid_[i]=func(grid_[i]);
554 4500 : for(unsigned j=0; j<dimension_; ++j) der_[i][j]=funcder(der_[i][j]);
555 : }
556 : } else {
557 0 : for(index_t i=0; i<grid_.size(); ++i) grid_[i]=func(grid_[i]);
558 : }
559 1 : }
560 :
561 970 : void Grid::writeHeader(OFile& ofile) {
562 3486 : for(unsigned i=0; i<dimension_; ++i) {
563 3774 : ofile.addConstantField("min_" + argnames[i]);
564 2516 : ofile.addConstantField("max_" + argnames[i]);
565 2516 : ofile.addConstantField("nbins_" + argnames[i]);
566 2516 : ofile.addConstantField("periodic_" + argnames[i]);
567 : }
568 970 : }
569 :
570 970 : void Grid::writeToFile(OFile& ofile) {
571 970 : vector<double> xx(dimension_);
572 970 : vector<double> der(dimension_);
573 : double f;
574 970 : writeHeader(ofile);
575 4802926 : for(index_t i=0; i<getSize(); ++i) {
576 4801956 : xx=getPoint(i);
577 2400978 : if(usederiv_) {f=getValueAndDerivatives(i,der);}
578 1863864 : else {f=getValue(i);}
579 7080881 : if(i>0 && dimension_>1 && getIndices(i)[dimension_-2]==0) ofile.printf("\n");
580 11761360 : for(unsigned j=0; j<dimension_; ++j) {
581 14040573 : ofile.printField("min_" + argnames[j], str_min_[j] );
582 9360382 : ofile.printField("max_" + argnames[j], str_max_[j] );
583 14040573 : ofile.printField("nbins_" + argnames[j], static_cast<int>(nbin_[j]) );
584 10761635 : if( pbc_[j] ) ofile.printField("periodic_" + argnames[j], "true" );
585 10612172 : else ofile.printField("periodic_" + argnames[j], "false" );
586 : }
587 21121742 : for(unsigned j=0; j<dimension_; ++j) { ofile.fmtField(" "+fmt_); ofile.printField(argnames[j],xx[j]); }
588 4801956 : ofile.fmtField(" "+fmt_); ofile.printField(funcname,f);
589 6617018 : if(usederiv_) for(unsigned j=0; j<dimension_; ++j) { ofile.fmtField(" "+fmt_); ofile.printField("der_" + argnames[j],der[j]); }
590 2400978 : ofile.printField();
591 : }
592 970 : }
593 :
594 0 : void Grid::writeCubeFile(OFile& ofile, const double& lunit) {
595 0 : plumed_assert( dimension_==3 );
596 0 : ofile.printf("PLUMED CUBE FILE\n");
597 0 : ofile.printf("OUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z\n");
598 : // Number of atoms followed by position of origin (origin set so that center of grid is in center of cell)
599 0 : ofile.printf("%d %f %f %f\n",1,-0.5*lunit*(max_[0]-min_[0]),-0.5*lunit*(max_[1]-min_[1]),-0.5*lunit*(max_[2]-min_[2]));
600 0 : ofile.printf("%u %f %f %f\n",nbin_[0],lunit*dx_[0],0.0,0.0); // Number of bins in each direction followed by
601 0 : ofile.printf("%u %f %f %f\n",nbin_[1],0.0,lunit*dx_[1],0.0); // shape of voxel
602 0 : ofile.printf("%u %f %f %f\n",nbin_[2],0.0,0.0,lunit*dx_[2]);
603 0 : ofile.printf("%d %f %f %f\n",1,0.0,0.0,0.0); // Fake atom otherwise VMD doesn't work
604 0 : std::vector<unsigned> pp(3);
605 0 : for(pp[0]=0; pp[0]<nbin_[0]; ++pp[0]) {
606 0 : for(pp[1]=0; pp[1]<nbin_[1]; ++pp[1]) {
607 0 : for(pp[2]=0; pp[2]<nbin_[2]; ++pp[2]) {
608 0 : ofile.printf("%f ",getValue(pp) );
609 0 : if(pp[2]%6==5) ofile.printf("\n");
610 : }
611 0 : ofile.printf("\n");
612 : }
613 : }
614 0 : }
615 :
616 18 : Grid* Grid::create(const std::string& funcl, const std::vector<Value*> & args, IFile& ifile,
617 : const vector<std::string> & gmin,const vector<std::string> & gmax,
618 : const vector<unsigned> & nbin,bool dosparse, bool dospline, bool doder) {
619 18 : Grid* grid=Grid::create(funcl,args,ifile,dosparse,dospline,doder);
620 18 : std::vector<unsigned> cbin( grid->getNbin() );
621 36 : std::vector<std::string> cmin( grid->getMin() ), cmax( grid->getMax() );
622 117 : for(unsigned i=0; i<args.size(); ++i) {
623 27 : plumed_massert( cmin[i]==gmin[i], "mismatched grid min" );
624 27 : plumed_massert( cmax[i]==gmax[i], "mismatched grid max" );
625 27 : if( args[i]->isPeriodic() ) {
626 16 : plumed_massert( cbin[i]==nbin[i], "mismatched grid nbins" );
627 : } else {
628 38 : plumed_massert( (cbin[i]-1)==nbin[i], "mismatched grid nbins");
629 : }
630 : }
631 18 : return grid;
632 : }
633 :
634 38 : Grid* Grid::create(const std::string& funcl, const std::vector<Value*> & args, IFile& ifile, bool dosparse, bool dospline, bool doder)
635 : {
636 : Grid* grid=NULL;
637 38 : unsigned nvar=args.size(); bool hasder=false; std::string pstring;
638 38 : std::vector<int> gbin1(nvar); std::vector<unsigned> gbin(nvar);
639 76 : std::vector<std::string> labels(nvar),gmin(nvar),gmax(nvar);
640 76 : std::vector<std::string> fieldnames; ifile.scanFieldList( fieldnames );
641 : // Retrieve names for fields
642 284 : for(unsigned i=0; i<args.size(); ++i) labels[i]=args[i]->getName();
643 : // And read the stuff from the header
644 38 : plumed_massert( ifile.FieldExist( funcl ), "no column labelled " + funcl + " in in grid input");
645 232 : for(unsigned i=0; i<args.size(); ++i) {
646 104 : ifile.scanField( "min_" + labels[i], gmin[i]);
647 104 : ifile.scanField( "max_" + labels[i], gmax[i]);
648 104 : ifile.scanField( "periodic_" + labels[i], pstring );
649 104 : ifile.scanField( "nbins_" + labels[i], gbin1[i]);
650 52 : plumed_assert( gbin1[i]>0 );
651 52 : if( args[i]->isPeriodic() ) {
652 18 : plumed_massert( pstring=="true", "input value is periodic but grid is not");
653 : std::string pmin, pmax;
654 54 : args[i]->getDomain( pmin, pmax ); gbin[i]=gbin1[i];
655 36 : if( pmin!=gmin[i] || pmax!=gmax[i] ) plumed_merror("mismatch between grid boundaries and periods of values");
656 : } else {
657 68 : gbin[i]=gbin1[i]-1; // Note header in grid file indicates one more bin that there should be when data is not periodic
658 34 : plumed_massert( pstring=="false", "input value is not periodic but grid is");
659 : }
660 156 : hasder=ifile.FieldExist( "der_" + args[i]->getName() );
661 52 : if( doder && !hasder ) plumed_merror("missing derivatives from grid file");
662 146 : for(unsigned j=0; j<fieldnames.size(); ++j) {
663 160 : for(unsigned k=i+1; k<args.size(); ++k) {
664 14 : if( fieldnames[j]==labels[k] ) plumed_merror("arguments in input are not in same order as in grid file");
665 : }
666 66 : if( fieldnames[j]==labels[i] ) break;
667 : }
668 : }
669 :
670 38 : if(!dosparse) {grid=new Grid(funcl,args,gmin,gmax,gbin,dospline,doder);}
671 0 : else {grid=new SparseGrid(funcl,args,gmin,gmax,gbin,dospline,doder);}
672 :
673 38 : vector<double> xx(nvar),dder(nvar);
674 38 : vector<double> dx=grid->getDx();
675 : double f,x;
676 93879 : while( ifile.scanField(funcl,f) ) {
677 464601 : for(unsigned i=0; i<nvar; ++i) {
678 741520 : ifile.scanField(labels[i],x); xx[i]=x+dx[i]/2.0;
679 370760 : ifile.scanField( "min_" + labels[i], gmin[i]);
680 370760 : ifile.scanField( "max_" + labels[i], gmax[i]);
681 370760 : ifile.scanField( "nbins_" + labels[i], gbin1[i]);
682 370760 : ifile.scanField( "periodic_" + labels[i], pstring );
683 : }
684 399742 : if(hasder) { for(unsigned i=0; i<nvar; ++i) { ifile.scanField( "der_" + args[i]->getName(), dder[i] ); } }
685 93841 : index_t index=grid->getIndex(xx);
686 93841 : if(doder) {grid->setValueAndDerivatives(index,f,dder);}
687 42609 : else {grid->setValue(index,f);}
688 93841 : ifile.scanField();
689 : }
690 38 : return grid;
691 : }
692 :
693 : // Sparse version of grid with map
694 0 : void SparseGrid::clear() {
695 : map_.clear();
696 0 : }
697 :
698 0 : Grid::index_t SparseGrid::getSize() const {
699 0 : return map_.size();
700 : }
701 :
702 0 : Grid::index_t SparseGrid::getMaxSize() const {
703 0 : return maxsize_;
704 : }
705 :
706 0 : double Grid::getDifferenceFromContour( const std::vector<double>& x, std::vector<double>& der ) const {
707 0 : return getValueAndDerivatives( x, der ) - contour_location;
708 : }
709 :
710 0 : void Grid::findSetOfPointsOnContour(const double& target, const std::vector<bool>& nosearch,
711 : unsigned& npoints, std::vector<std::vector<double> >& points ) {
712 : // Set contour location for function
713 0 : contour_location=target;
714 : // Resize points to maximum possible value
715 0 : points.resize( dimension_*maxsize_ );
716 :
717 : // Two points for search
718 0 : std::vector<unsigned> ind(dimension_);
719 0 : std::vector<double> direction( dimension_, 0 );
720 :
721 : // Run over whole grid
722 0 : npoints=0; RootFindingBase<Grid> mymin( this );
723 0 : for(unsigned i=0; i<maxsize_; ++i) {
724 0 : for(unsigned j=0; j<dimension_; ++j) ind[j]=getIndices(i)[j];
725 :
726 : // Get the value of a point on the grid
727 0 : double val1=getValue(i) - target;
728 :
729 : // Now search for contour in each direction
730 : bool edge=false;
731 0 : for(unsigned j=0; j<dimension_; ++j) {
732 0 : if( nosearch[j] ) continue ;
733 : // Make sure we don't search at the edge of the grid
734 0 : if( !pbc_[j] && (ind[j]+1)==nbin_[j] ) continue;
735 0 : else if( (ind[j]+1)==nbin_[j] ) { edge=true; ind[j]=0; }
736 0 : else ind[j]+=1;
737 0 : double val2=getValue(ind) - target;
738 0 : if( val1*val2<0 ) {
739 : // Use initial point location as first guess for search
740 0 : points[npoints].resize(dimension_); for(unsigned k=0; k<dimension_; ++k) points[npoints][k]=getPoint(i)[k];
741 : // Setup direction vector
742 0 : direction[j]=0.999999999*dx_[j];
743 : // And do proper search for contour point
744 0 : mymin.linesearch( direction, points[npoints], &Grid::getDifferenceFromContour );
745 0 : direction[j]=0.0; npoints++;
746 : }
747 0 : if( pbc_[j] && edge ) { edge=false; ind[j]=nbin_[j]-1; }
748 0 : else ind[j]-=1;
749 : }
750 : }
751 0 : }
752 :
753 0 : double SparseGrid::getValue(index_t index)const {
754 0 : plumed_assert(index<maxsize_);
755 : double value=0.0;
756 : const auto it=map_.find(index);
757 0 : if(it!=map_.end()) value=it->second;
758 0 : return value;
759 : }
760 :
761 440 : double SparseGrid::getValueAndDerivatives
762 : (index_t index, vector<double>& der)const {
763 880 : plumed_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
764 : double value=0.0;
765 2080 : for(unsigned int i=0; i<dimension_; ++i) der[i]=0.0;
766 : const auto it=map_.find(index);
767 440 : if(it!=map_.end()) value=it->second;
768 : const auto itder=der_.find(index);
769 440 : if(itder!=der_.end()) der=itder->second;
770 440 : return value;
771 : }
772 :
773 0 : void SparseGrid::setValue(index_t index, double value) {
774 0 : plumed_assert(index<maxsize_ && !usederiv_);
775 0 : map_[index]=value;
776 0 : }
777 :
778 0 : void SparseGrid::setValueAndDerivatives
779 : (index_t index, double value, vector<double>& der) {
780 0 : plumed_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
781 0 : map_[index]=value;
782 0 : der_[index]=der;
783 0 : }
784 :
785 0 : void SparseGrid::addValue(index_t index, double value) {
786 0 : plumed_assert(index<maxsize_ && !usederiv_);
787 0 : map_[index]+=value;
788 0 : }
789 :
790 19999 : void SparseGrid::addValueAndDerivatives
791 : (index_t index, double value, vector<double>& der) {
792 39998 : plumed_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
793 19999 : map_[index]+=value;
794 19999 : der_[index].resize(dimension_);
795 139048 : for(unsigned int i=0; i<dimension_; ++i) der_[index][i]+=der[i];
796 19999 : }
797 :
798 0 : void SparseGrid::writeToFile(OFile& ofile) {
799 0 : vector<double> xx(dimension_);
800 0 : vector<double> der(dimension_);
801 : double f;
802 0 : writeHeader(ofile);
803 0 : ofile.fmtField(" "+fmt_);
804 0 : for(const auto & it : map_) {
805 0 : index_t i=it.first;
806 0 : xx=getPoint(i);
807 0 : if(usederiv_) {f=getValueAndDerivatives(i,der);}
808 0 : else {f=getValue(i);}
809 0 : if(i>0 && dimension_>1 && getIndices(i)[dimension_-2]==0) ofile.printf("\n");
810 0 : for(unsigned j=0; j<dimension_; ++j) {
811 0 : ofile.printField("min_" + argnames[j], str_min_[j] );
812 0 : ofile.printField("max_" + argnames[j], str_max_[j] );
813 0 : ofile.printField("nbins_" + argnames[j], static_cast<int>(nbin_[j]) );
814 0 : if( pbc_[j] ) ofile.printField("periodic_" + argnames[j], "true" );
815 0 : else ofile.printField("periodic_" + argnames[j], "false" );
816 : }
817 0 : for(unsigned j=0; j<dimension_; ++j) ofile.printField(argnames[j],xx[j]);
818 0 : ofile.printField(funcname, f);
819 0 : if(usederiv_) { for(unsigned j=0; j<dimension_; ++j) ofile.printField("der_" + argnames[j],der[j]); }
820 0 : ofile.printField();
821 : }
822 0 : }
823 :
824 0 : double SparseGrid::getMinValue() const {
825 : double minval;
826 : minval=0.0;
827 0 : for(auto const & i : map_) {
828 0 : if(i.second<minval) minval=i.second;
829 : }
830 0 : return minval;
831 : }
832 :
833 9 : double SparseGrid::getMaxValue() const {
834 : double maxval;
835 : maxval=0.0;
836 590 : for(auto const & i : map_) {
837 581 : if(i.second>maxval) maxval=i.second;
838 : }
839 9 : return maxval;
840 : }
841 :
842 :
843 688228 : void Grid::projectOnLowDimension(double &val, std::vector<int> &vHigh, WeightBase * ptr2obj ) {
844 : unsigned i=0;
845 5475683 : for(i=0; i<vHigh.size(); i++) {
846 1373115 : if(vHigh[i]<0) { // this bin needs to be integrated out
847 : // parallelize here???
848 2746206 : for(unsigned j=0; j<(getNbin())[i]; j++) {
849 681522 : vHigh[i]=int(j);
850 681522 : projectOnLowDimension(val,vHigh,ptr2obj); // recursive function: this is the core of the mechanism
851 681522 : vHigh[i]=-1;
852 : }
853 : return; //
854 : }
855 : }
856 : // when there are no more bin to dig in then retrieve the value
857 681522 : if(i==vHigh.size()) {
858 : //std::cerr<<"POINT: ";
859 : //for(unsigned j=0;j<vHigh.size();j++){
860 : // std::cerr<<vHigh[j]<<" ";
861 : //}
862 681522 : std::vector<unsigned> vv(vHigh.size());
863 6815220 : for(unsigned j=0; j<vHigh.size(); j++)vv[j]=unsigned(vHigh[j]);
864 : //
865 : // this is the real assignment !!!!! (hack this to have bias or other stuff)
866 : //
867 :
868 : // this case: produce fes
869 : //val+=exp(beta*getValue(vv)) ;
870 681522 : double myv=getValue(vv);
871 681522 : val=ptr2obj->projectInnerLoop(val,myv) ;
872 : // to be added: bias (same as before without negative sign)
873 : //std::cerr<<" VAL: "<<val <<endl;
874 : }
875 : }
876 :
877 66 : Grid Grid::project(const std::vector<std::string> & proj, WeightBase *ptr2obj ) {
878 : // find extrema only for the projection
879 66 : vector<string> smallMin,smallMax;
880 : vector<unsigned> smallBin;
881 : vector<unsigned> dimMapping;
882 : vector<bool> smallIsPeriodic;
883 66 : vector<string> smallName;
884 :
885 : // check if the two key methods are there
886 : WeightBase* pp = dynamic_cast<WeightBase*>(ptr2obj);
887 66 : if (!pp)plumed_merror("This WeightBase is not complete: you need a projectInnerLoop and projectOuterLoop ");
888 :
889 330 : for(unsigned j=0; j<proj.size(); j++) {
890 196 : for(unsigned i=0; i<getArgNames().size(); i++) {
891 196 : if(proj[j]==getArgNames()[i]) {
892 : unsigned offset;
893 : // note that at sizetime the non periodic dimension get a bin more
894 132 : if(getIsPeriodic()[i]) {offset=0;} else {offset=1;}
895 132 : smallMax.push_back(getMax()[i]);
896 132 : smallMin.push_back(getMin()[i]);
897 264 : smallBin.push_back(getNbin()[i]-offset);
898 198 : smallIsPeriodic.push_back(getIsPeriodic()[i]);
899 66 : dimMapping.push_back(i);
900 132 : smallName.push_back(getArgNames()[i]);
901 66 : break;
902 : }
903 : }
904 : }
905 132 : Grid smallgrid("projection",smallName,smallMin,smallMax,smallBin,false,false,true,smallIsPeriodic,smallMin,smallMax);
906 : // check that the two grids are commensurate
907 330 : for(unsigned i=0; i<dimMapping.size(); i++) {
908 198 : plumed_massert( (smallgrid.getMax())[i] == (getMax())[dimMapping[i]], "the two input grids are not compatible in max" );
909 198 : plumed_massert( (smallgrid.getMin())[i] == (getMin())[dimMapping[i]], "the two input grids are not compatible in min" );
910 396 : plumed_massert( (smallgrid.getNbin())[i]== (getNbin())[dimMapping[i]], "the two input grids are not compatible in bin" );
911 : }
912 : vector<unsigned> toBeIntegrated;
913 396 : for(unsigned i=0; i<getArgNames().size(); i++) {
914 : bool doappend=true;
915 462 : for(unsigned j=0; j<dimMapping.size(); j++) {
916 132 : if(dimMapping[j]==i) {doappend=false; break;}
917 : }
918 132 : if(doappend)toBeIntegrated.push_back(i);
919 : }
920 : //for(unsigned i=0;i<dimMapping.size();i++ ){
921 : // cerr<<"Dimension to preserve "<<dimMapping[i]<<endl;
922 : //}
923 : //for(unsigned i=0;i<toBeIntegrated.size();i++ ){
924 : // cerr<<"Dimension to integrate "<<toBeIntegrated[i]<<endl;
925 : //}
926 :
927 : // loop over all the points in the Grid, find the corresponding fixed index, rotate over all the other ones
928 13478 : for(unsigned i=0; i<smallgrid.getSize(); i++) {
929 : std::vector<unsigned> v;
930 13412 : v=smallgrid.getIndices(i);
931 13412 : std::vector<int> vHigh((getArgNames()).size(),-1);
932 46942 : for(unsigned j=0; j<dimMapping.size(); j++)vHigh[dimMapping[j]]=int(v[j]);
933 : // the vector vhigh now contains at the beginning the index of the low dimension and -1 for the values that need to be integrated
934 6706 : double initval=0.;
935 6706 : projectOnLowDimension(initval,vHigh, ptr2obj);
936 6706 : smallgrid.setValue(i,initval);
937 : }
938 : // reset to zero just for biasing (this option can be evtl enabled in a future...)
939 : //double vmin;vmin=-smallgrid.getMinValue()+1;
940 13478 : for(unsigned i=0; i<smallgrid.getSize(); i++) {
941 : // //if(dynamic_cast<BiasWeight*>(ptr2obj)){
942 : // // smallgrid.addValue(i,vmin);// go to 1
943 : // //}
944 6706 : double vv=smallgrid.getValue(i);
945 6706 : smallgrid.setValue(i,ptr2obj->projectOuterLoop(vv));
946 : // //if(dynamic_cast<BiasWeight*>(ptr2obj)){
947 : // // smallgrid.addValue(i,-vmin);// bring back to the value
948 : // //}
949 : }
950 :
951 66 : return smallgrid;
952 : }
953 :
954 0 : double Grid::integrate( std::vector<unsigned>& npoints ) {
955 0 : plumed_dbg_assert( npoints.size()==dimension_ ); plumed_assert( dospline_ );
956 :
957 : unsigned ntotgrid=1; double box_vol=1.0;
958 0 : std::vector<double> ispacing( npoints.size() );
959 0 : for(unsigned j=0; j<dimension_; ++j) {
960 0 : if( !pbc_[j] ) {
961 0 : ispacing[j] = ( max_[j] - dx_[j] - min_[j] ) / static_cast<double>( npoints[j] );
962 0 : npoints[j]+=1;
963 : } else {
964 0 : ispacing[j] = ( max_[j] - min_[j] ) / static_cast<double>( npoints[j] );
965 : }
966 0 : ntotgrid*=npoints[j]; box_vol*=ispacing[j];
967 : }
968 :
969 0 : std::vector<double> vals( dimension_ );
970 0 : std::vector<unsigned> t_index( dimension_ ); double integral=0.0;
971 0 : for(unsigned i=0; i<ntotgrid; ++i) {
972 0 : t_index[0]=(i%npoints[0]);
973 : unsigned kk=i;
974 0 : for(unsigned j=1; j<dimension_-1; ++j) { kk=(kk-t_index[j-1])/npoints[i-1]; t_index[j]=(kk%npoints[i]); }
975 0 : if( dimension_>=2 ) t_index[dimension_-1]=((kk-t_index[dimension_-1])/npoints[dimension_-2]);
976 :
977 0 : for(unsigned j=0; j<dimension_; ++j) vals[j]=min_[j] + t_index[j]*ispacing[j];
978 :
979 0 : integral += getValue( vals );
980 : }
981 :
982 0 : return box_vol*integral;
983 : }
984 :
985 0 : void Grid::mpiSumValuesAndDerivatives( Communicator& comm ) {
986 0 : comm.Sum( grid_ ); for(unsigned i=0; i<der_.size(); ++i) comm.Sum( der_[i] );
987 0 : }
988 :
989 :
990 59582 : bool indexed_lt(pair<Grid::index_t, double> const &x, pair<Grid::index_t, double> const &y) {
991 59582 : return x.second < y.second;
992 : }
993 :
994 273 : double Grid::findMaximalPathMinimum(const std::vector<double> &source, const std::vector<double> &sink) {
995 : plumed_dbg_assert(source.size() == dimension_);
996 : plumed_dbg_assert(sink.size() == dimension_);
997 : // Start and end indices
998 273 : index_t source_idx = getIndex(source);
999 273 : index_t sink_idx = getIndex(sink);
1000 : // Path cost
1001 : double maximal_minimum = 0;
1002 : // In one dimension, path searching is very easy--either go one way if it's not periodic,
1003 : // or go both ways if it is periodic. There's no reason to pay the cost of Dijkstra.
1004 273 : if (dimension_ == 1) {
1005 : // Do a search from the grid source to grid sink that does not
1006 : // cross the grid boundary.
1007 147 : double curr_min_bias = getValue(source_idx);
1008 : // Either search from a high source to a low sink.
1009 147 : if (source_idx > sink_idx) {
1010 1323 : for (index_t i = source_idx; i >= sink_idx; i--) {
1011 588 : if (curr_min_bias == 0.0) {
1012 : break;
1013 : }
1014 588 : curr_min_bias = fmin(curr_min_bias, getValue(i));
1015 : }
1016 : // Or search from a low source to a high sink.
1017 0 : } else if (source_idx < sink_idx) {
1018 0 : for (index_t i = source_idx; i <= sink_idx; i++) {
1019 0 : if (curr_min_bias == 0.0) {
1020 : break;
1021 : }
1022 0 : curr_min_bias = fmin(curr_min_bias, getValue(i));
1023 : }
1024 : }
1025 : maximal_minimum = curr_min_bias;
1026 : // If the grid is periodic, also do the search that crosses
1027 : // the grid boundary.
1028 147 : if (pbc_[0]) {
1029 42 : double curr_min_bias = getValue(source_idx);
1030 : // Either go from a high source to the upper boundary and
1031 : // then from the bottom boundary to the sink
1032 42 : if (source_idx > sink_idx) {
1033 378 : for (index_t i = source_idx; i < maxsize_; i++) {
1034 168 : if (curr_min_bias == 0.0) {
1035 : break;
1036 : }
1037 168 : curr_min_bias = fmin(curr_min_bias, getValue(i));
1038 : }
1039 378 : for (index_t i = 0; i <= sink_idx; i++) {
1040 168 : if (curr_min_bias == 0.0) {
1041 : break;
1042 : }
1043 168 : curr_min_bias = fmin(curr_min_bias, getValue(i));
1044 : }
1045 : // Or go from a low source to the bottom boundary and
1046 : // then from the high boundary to the sink
1047 0 : } else if (source_idx < sink_idx) {
1048 0 : for (index_t i = source_idx; i > 0; i--) {
1049 0 : if (curr_min_bias == 0.0) {
1050 : break;
1051 : }
1052 0 : curr_min_bias = fmin(curr_min_bias, getValue(i));
1053 : }
1054 0 : curr_min_bias = fmin(curr_min_bias, getValue(0));
1055 0 : for (index_t i = maxsize_ - 1; i <= sink_idx; i--) {
1056 0 : if (curr_min_bias == 0.0) {
1057 : break;
1058 : }
1059 0 : curr_min_bias = fmin(curr_min_bias, getValue(i));
1060 : }
1061 : }
1062 : // If the boundary crossing paths was more biased, it's
1063 : // minimal bias replaces the non-boundary-crossing path's
1064 : // minimum.
1065 42 : maximal_minimum = fmax(maximal_minimum, curr_min_bias);
1066 : }
1067 : // The one dimensional path search is complete.
1068 : return maximal_minimum;
1069 : // In two or more dimensions, path searching isn't trivial and we really
1070 : // do need to use a path search algorithm. Dijkstra is the simplest decent
1071 : // one. Using it we've never found the path search to be performance
1072 : // limiting in any solvated biomolecule test system, but faster options are
1073 : // easy to imagine if they become necessary. NB-In this case, we're actually
1074 : // using a greedy variant of Dijkstra's algorithm where the first possible
1075 : // path to a point always controls the path cost to that point. The structure
1076 : // of the cost function in this case guarantees that the calculated costs will
1077 : // be correct using this variant even though fine details of the paths may not
1078 : // match a normal Dijkstra search.
1079 126 : } else if (dimension_ > 1) {
1080 : // Prepare calculation temporaries for Dijkstra's algorithm.
1081 : // Minimal path costs from source to a given grid point
1082 126 : vector<double> mins_from_source = vector<double>(maxsize_, -1.0);
1083 : // Heap for tracking available steps, steps are recorded as std::pairs of
1084 : // an index and a value.
1085 : vector< pair<index_t, double> > next_steps;
1086 : pair<index_t, double> curr_indexed_val;
1087 : make_heap(next_steps.begin(), next_steps.end(), indexed_lt);
1088 : // The search begins at the source index.
1089 252 : next_steps.push_back(pair<index_t, double>(source_idx, getValue(source_idx)));
1090 126 : push_heap(next_steps.begin(), next_steps.end(), indexed_lt);
1091 : // At first no points have been examined and the optimal path has not been found.
1092 : index_t n_examined = 0;
1093 : bool path_not_found = true;
1094 : // Until a path is found,
1095 9576 : while (path_not_found) {
1096 : // Examine the grid point currently most accessible from
1097 : // the set of all previously explored grid points by popping
1098 : // it from the top of the heap.
1099 : pop_heap(next_steps.begin(), next_steps.end(), indexed_lt);
1100 : curr_indexed_val = next_steps.back();
1101 : next_steps.pop_back();
1102 : n_examined++;
1103 : // Check if this point is the sink point, and if so
1104 : // finish the loop.
1105 9702 : if (curr_indexed_val.first == sink_idx) {
1106 : path_not_found = false;
1107 : maximal_minimum = curr_indexed_val.second;
1108 126 : break;
1109 : // Check if this point has reached the worst possible
1110 : // value, and if so stop looking for paths.
1111 9576 : } else if (curr_indexed_val.second == 0.0) {
1112 : maximal_minimum = 0.0;
1113 : break;
1114 : }
1115 : // If the search is not over, add this grid point's neighbors to the
1116 : // possible next points to search for the sink.
1117 9576 : vector<index_t> neighs = getNearestNeighbors(curr_indexed_val.first);
1118 129162 : for (unsigned k = 0; k < neighs.size(); k++) {
1119 36670 : index_t i = neighs[k];
1120 : // If the neighbor has not already been added to the list of possible next steps,
1121 36670 : if (mins_from_source[i] == -1.0) {
1122 : // Set the cost to reach it via a path through the current point being examined.
1123 24568 : mins_from_source[i] = fmin(curr_indexed_val.second, getValue(i));
1124 : // Add the neighboring point to the heap of potential next steps.
1125 12284 : next_steps.push_back(pair<index_t, double>(i, mins_from_source[i]));
1126 12284 : push_heap(next_steps.begin(), next_steps.end(), indexed_lt);
1127 : }
1128 : }
1129 : // Move on to the next best looking step along any of the paths
1130 : // growing from the source.
1131 : }
1132 : // The multidimensional path search is now complete.
1133 : return maximal_minimum;
1134 : }
1135 : return 0.0;
1136 : }
1137 :
1138 4839 : }
|
