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