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