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