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