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