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