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 <limits>
35 : #include <cstdlib>
36 : #include <cmath>
37 : #include <iostream>
38 : #include <sstream>
39 : #include <cstdio>
40 : #include <cfloat>
41 : #include <array>
42 :
43 : namespace PLMD {
44 :
45 : constexpr std::size_t GridBase::maxdim;
46 :
47 : template<unsigned dimension_>
48 3242 : class Accelerator final:
49 : public Grid::AcceleratorBase {
50 :
51 : public:
52 :
53 1654 : unsigned getDimension() const override {
54 1654 : return dimension_;
55 : }
56 :
57 61864 : 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 {
58 : plumed_dbg_assert(indices_size==dimension_ && nneigh.size()==dimension_);
59 :
60 61864 : std::vector<Grid::index_t> neighbors;
61 : std::array<unsigned,dimension_> small_bin;
62 : std::array<unsigned,dimension_> small_indices;
63 : std::array<unsigned,dimension_> tmp_indices;
64 :
65 : unsigned small_nbin=1;
66 184870 : for(unsigned j=0; j<dimension_; ++j) {
67 123006 : small_bin[j]=(2*nneigh[j]+1);
68 : // along a periodic direction a neighborhood larger than the grid would visit
69 : // some points more than once: visit each point of the grid exactly once instead
70 123006 : if(pbc_[j] && small_bin[j]>nbin_[j]) {
71 0 : small_bin[j]=nbin_[j];
72 : }
73 123006 : small_nbin*=small_bin[j];
74 : }
75 61864 : neighbors.reserve(small_nbin);
76 :
77 184870 : for(unsigned j=0; j<dimension_; j++) {
78 123006 : small_indices[j]=0;
79 : }
80 :
81 1373282 : for(unsigned index=0; index<small_nbin; ++index) {
82 : unsigned ll=0;
83 3768576 : for(unsigned i=0; i<dimension_; ++i) {
84 2457158 : int i0=small_indices[i]-nneigh[i]+indices[i];
85 2457158 : if(!pbc_[i] && i0<0) {
86 46880 : continue;
87 : }
88 2410278 : if(!pbc_[i] && i0>=static_cast<int>(nbin_[i])) {
89 35216 : continue;
90 : }
91 2375062 : if( pbc_[i] ) {
92 2116244 : const int n=nbin_[i];
93 2116244 : i0=((i0%n)+n)%n;
94 : }
95 2375062 : tmp_indices[ll]=static_cast<unsigned>(i0);
96 2375062 : ll++;
97 : }
98 1311418 : if(ll==dimension_) {
99 1229322 : neighbors.push_back(getIndex(grid,nbin_,&tmp_indices[0],dimension_));
100 : }
101 :
102 1311418 : small_indices[0]++;
103 2291480 : for(unsigned j=0; j<dimension_-1; j++)
104 1145740 : if(small_indices[j]==small_bin[j]) {
105 162562 : small_indices[j]=0;
106 162562 : small_indices[j+1]++;
107 : }
108 : }
109 61864 : return neighbors;
110 : }
111 :
112 : // we are flattening arrays using a column-major order
113 7293498 : GridBase::index_t getIndex(const GridBase& grid, const std::vector<unsigned> & nbin_,const unsigned* indices,std::size_t indices_size) const override {
114 19634249 : for(unsigned int i=0; i<dimension_; i++)
115 12340751 : if(indices[i]>=nbin_[i]) {
116 0 : plumed_error() << "Looking for a value outside the grid along the " << i << " dimension (arg name: "<<grid.getArgNames()[i]<<")";
117 : }
118 7293498 : auto index=indices[dimension_-1];
119 9658156 : for(unsigned int i=dimension_-1; i>0; --i) {
120 5047253 : index=index*nbin_[i-1]+indices[i-1];
121 : }
122 7293498 : return index;
123 : }
124 :
125 33313934 : 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 {
126 113361970 : for(unsigned int i=0; i<dimension_; ++i) {
127 80048036 : point[i]=min_[i]+(double)(indices[i])*dx_[i];
128 : }
129 33313934 : }
130 :
131 1203510 : 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 {
132 : plumed_dbg_assert(x.size()==dimension_);
133 : plumed_dbg_assert(rindex_size==dimension_);
134 3566946 : for(unsigned int i=0; i<dimension_; ++i) {
135 2363436 : const double f=std::floor((x[i]-min_[i])/dx_[i]);
136 : // converting a negative (or too large) double to unsigned is undefined behavior:
137 : // on x86 it gave a huge index (and getIndex raised an error), on arm64 it gives 0
138 : // (i.e. the value of the first bin). Map it explicitly to an invalid index.
139 2363436 : rindex_data[i] = (f>=0.0 && f<double(std::numeric_limits<unsigned>::max())) ?
140 : static_cast<unsigned>(f) : std::numeric_limits<unsigned>::max();
141 : }
142 1203510 : }
143 :
144 :
145 : // we are flattening arrays using a column-major order
146 38013328 : void getIndices(const std::vector<unsigned> & nbin_, GridBase::index_t index, unsigned* indices, std::size_t indices_size) const override {
147 38013328 : plumed_assert(indices_size==dimension_)<<indices_size;
148 91039369 : for(unsigned int i=0; i<dimension_-1; ++i) {
149 53994039 : indices[i]=index%nbin_[i];
150 53994039 : index/=nbin_[i];
151 : }
152 38013328 : indices[dimension_-1]=index;
153 38013328 : }
154 :
155 : };
156 :
157 3242 : std::unique_ptr<Grid::AcceleratorBase> Grid::AcceleratorBase::create(unsigned dim) {
158 : // I do this by enumeration.
159 : // Maybe can be simplified using the preprocessor, but I am not sure it's worth it
160 : if(dim==1) {
161 2573 : return std::make_unique<Accelerator<1>>();
162 : }
163 : if(dim==2) {
164 661 : return std::make_unique<Accelerator<2>>();
165 : }
166 : if(dim==3) {
167 8 : return std::make_unique<Accelerator<3>>();
168 : }
169 : if(dim==4) {
170 0 : return std::make_unique<Accelerator<4>>();
171 : }
172 : if(dim==5) {
173 0 : return std::make_unique<Accelerator<5>>();
174 : }
175 : if(dim==6) {
176 0 : return std::make_unique<Accelerator<6>>();
177 : }
178 : if(dim==7) {
179 0 : return std::make_unique<Accelerator<7>>();
180 : }
181 : if(dim==8) {
182 0 : return std::make_unique<Accelerator<8>>();
183 : }
184 : if(dim==9) {
185 0 : return std::make_unique<Accelerator<9>>();
186 : }
187 : if(dim==10) {
188 0 : return std::make_unique<Accelerator<10>>();
189 : }
190 : if(dim==11) {
191 0 : return std::make_unique<Accelerator<11>>();
192 : }
193 : if(dim==12) {
194 0 : return std::make_unique<Accelerator<12>>();
195 : }
196 : if(dim==13) {
197 0 : return std::make_unique<Accelerator<13>>();
198 : }
199 : if(dim==14) {
200 0 : return std::make_unique<Accelerator<14>>();
201 : }
202 : if(dim==15) {
203 0 : return std::make_unique<Accelerator<15>>();
204 : }
205 : if(dim==16) {
206 0 : return std::make_unique<Accelerator<16>>();
207 : }
208 : // no need to go beyond this
209 0 : plumed_error();
210 : }
211 :
212 1460 : GridBase::GridBase(const std::string& funcl, const std::vector<Value*> & args, const std::vector<std::string> & gmin,
213 1460 : const std::vector<std::string> & gmax, const std::vector<unsigned> & nbin, bool dospline, bool usederiv) {
214 : // various checks
215 1460 : plumed_assert(args.size()<=maxdim) << "grid dim cannot exceed "<<maxdim;
216 1460 : plumed_massert(args.size()==gmin.size(),"grid min dimensions in input do not match number of arguments");
217 1460 : plumed_massert(args.size()==nbin.size(),"number of bins on input do not match number of arguments");
218 1460 : plumed_massert(args.size()==gmax.size(),"grid max dimensions in input do not match number of arguments");
219 1460 : unsigned dim=gmax.size();
220 1460 : std::vector<std::string> names;
221 1460 : std::vector<bool> isperiodic;
222 1460 : std::vector<std::string> pmin,pmax;
223 1460 : names.resize( dim );
224 1460 : isperiodic.resize( dim );
225 1460 : pmin.resize( dim );
226 1460 : pmax.resize( dim );
227 3248 : for(unsigned int i=0; i<dim; ++i) {
228 1788 : names[i]=args[i]->getName();
229 1788 : if( args[i]->isPeriodic() ) {
230 : isperiodic[i]=true;
231 595 : args[i]->getDomain( pmin[i], pmax[i] );
232 : } else {
233 : isperiodic[i]=false;
234 : pmin[i]="0.";
235 : pmax[i]="0.";
236 : }
237 : }
238 : // this is a value-independent initializator
239 1460 : Init(funcl,names,gmin,gmax,nbin,dospline,usederiv,isperiodic,pmin,pmax);
240 2920 : }
241 :
242 128 : GridBase::GridBase(const std::string& funcl, const std::vector<std::string> &names, const std::vector<std::string> & gmin,
243 : const std::vector<std::string> & gmax, const std::vector<unsigned> & nbin, bool dospline, bool usederiv,
244 128 : const std::vector<bool> &isperiodic, const std::vector<std::string> &pmin, const std::vector<std::string> &pmax ) {
245 : // this calls the initializator
246 128 : Init(funcl,names,gmin,gmax,nbin,dospline,usederiv,isperiodic,pmin,pmax);
247 128 : }
248 :
249 1588 : void GridBase::Init(const std::string& funcl, const std::vector<std::string> &names, const std::vector<std::string> & gmin,
250 : const std::vector<std::string> & gmax, const std::vector<unsigned> & nbin, bool dospline, bool usederiv,
251 : const std::vector<bool> &isperiodic, const std::vector<std::string> &pmin, const std::vector<std::string> &pmax ) {
252 1588 : fmt_="%14.9f";
253 : // various checks
254 1588 : plumed_assert(names.size()<=maxdim) << "grid size cannot exceed "<<maxdim;
255 1588 : plumed_massert(names.size()==gmin.size(),"grid dimensions in input do not match number of arguments");
256 1588 : plumed_massert(names.size()==nbin.size(),"grid dimensions in input do not match number of arguments");
257 1588 : plumed_massert(names.size()==gmax.size(),"grid dimensions in input do not match number of arguments");
258 1588 : dimension_=gmax.size();
259 1588 : accelerator=AcceleratorHandler(dimension_);
260 1588 : str_min_=gmin;
261 1588 : str_max_=gmax;
262 1588 : argnames.resize( dimension_ );
263 1588 : min_.resize( dimension_ );
264 1588 : max_.resize( dimension_ );
265 1588 : pbc_.resize( dimension_ );
266 3504 : for(unsigned int i=0; i<dimension_; ++i) {
267 1916 : argnames[i]=names[i];
268 1916 : if( isperiodic[i] ) {
269 : pbc_[i]=true;
270 : str_min_[i]=pmin[i];
271 : str_max_[i]=pmax[i];
272 : } else {
273 : pbc_[i]=false;
274 : }
275 1916 : Tools::convert(str_min_[i],min_[i]);
276 1916 : Tools::convert(str_max_[i],max_[i]);
277 1916 : funcname=funcl;
278 1916 : plumed_massert(max_[i]>min_[i],"maximum in grid must be larger than minimum");
279 1916 : plumed_massert(nbin[i]>0,"number of grid points must be greater than zero");
280 : }
281 1588 : nbin_=nbin;
282 1588 : dospline_=dospline;
283 1588 : usederiv_=usederiv;
284 1588 : if(dospline_) {
285 216 : plumed_assert(dospline_==usederiv_);
286 : }
287 1588 : maxsize_=1;
288 3504 : for(unsigned int i=0; i<dimension_; ++i) {
289 1916 : dx_.push_back( (max_[i]-min_[i])/static_cast<double>( nbin_[i] ) );
290 1916 : if( !pbc_[i] ) {
291 1257 : max_[i] += dx_[i];
292 1257 : nbin_[i] += 1;
293 : }
294 1916 : maxsize_*=nbin_[i];
295 : }
296 1588 : }
297 :
298 45997 : bool GridBase::isInside(const std::vector<double> & x) const {
299 : plumed_dbg_assert(x.size()==dimension_);
300 110272 : for(unsigned i=0; i<dimension_; ++i) {
301 64603 : if(!pbc_[i]) {
302 : // same arithmetic as getIndices (false also for a NaN)
303 30043 : const double k=std::floor((x[i]-min_[i])/dx_[i]);
304 30043 : if(!(k>=0.0 && k<static_cast<double>(nbin_[i]))) {
305 : return false;
306 : }
307 : }
308 : }
309 : return true;
310 : }
311 :
312 : namespace {
313 : /// bound of a resized grid as a string: 12 significant digits, so that the files do not depend on the last digits of
314 : /// the bounds, which can differ between platforms (e.g. with fused multiply-add); when a grid is read from a file
315 : /// its lattice is compared with a tolerance
316 266 : std::string boundString(double x) {
317 : char buffer[64];
318 : std::snprintf(buffer,sizeof(buffer),"%.12g",x);
319 266 : return buffer;
320 : }
321 : }
322 :
323 132 : void GridBase::resize(const std::vector<long long> & lower, const std::vector<long long> & upper) {
324 132 : plumed_assert(lower.size()==dimension_ && upper.size()==dimension_);
325 132 : const std::vector<unsigned> oldnbin(nbin_);
326 132 : std::vector<unsigned> newnbin(nbin_);
327 : double newsize=1.0;
328 318 : for(unsigned i=0; i<dimension_; ++i) {
329 186 : if(lower[i]!=0 || upper[i]!=0) {
330 133 : plumed_massert(!pbc_[i],"periodic directions of a grid cannot be resized");
331 133 : const long long n=static_cast<long long>(nbin_[i])+lower[i]+upper[i];
332 133 : plumed_massert(n>=2 && n<=static_cast<long long>(std::numeric_limits<unsigned>::max()),"invalid number of grid points after resizing");
333 133 : newnbin[i]=static_cast<unsigned>(n);
334 : }
335 186 : newsize*=newnbin[i];
336 : }
337 132 : plumed_massert(newsize<static_cast<double>(std::numeric_limits<index_t>::max()),"too many grid points after resizing");
338 318 : for(unsigned i=0; i<dimension_; ++i) {
339 186 : if(lower[i]!=0 || upper[i]!=0) {
340 133 : min_[i]-=static_cast<double>(lower[i])*dx_[i];
341 133 : max_[i]+=static_cast<double>(upper[i])*dx_[i];
342 133 : nbin_[i]=newnbin[i];
343 : // as in the input of a non-periodic grid, the maximum is the last point
344 266 : str_min_[i]=boundString(min_[i]);
345 266 : str_max_[i]=boundString(max_[i]-dx_[i]);
346 : }
347 : }
348 132 : maxsize_=1;
349 318 : for(unsigned i=0; i<dimension_; ++i) {
350 186 : maxsize_*=nbin_[i];
351 : }
352 132 : resizeData(oldnbin,lower);
353 132 : }
354 :
355 567 : void GridBase::forEachStoredPoint(const std::function<void(index_t,double)> & f) const {
356 567 : const index_t n=getSize();
357 9412884 : for(index_t i=0; i<n; ++i) {
358 9412317 : f(i,getValue(i));
359 : }
360 567 : }
361 :
362 436 : void SparseGrid::forEachStoredPoint(const std::function<void(index_t,double)> & f) const {
363 1095308 : for(const auto & it : map_) {
364 1094872 : f(it.first,it.second);
365 : }
366 436 : }
367 :
368 0 : void GridBase::resizeData(const std::vector<unsigned> & oldnbin, const std::vector<long long> & shift) {
369 0 : plumed_merror("this kind of grid cannot be resized");
370 : }
371 :
372 1360671 : std::vector<std::string> GridBase::getMin() const {
373 1360671 : return str_min_;
374 : }
375 :
376 559 : std::vector<std::string> GridBase::getMax() const {
377 559 : return str_max_;
378 : }
379 :
380 1330206 : std::vector<double> GridBase::getDx() const {
381 1330206 : return dx_;
382 : }
383 :
384 213050 : double GridBase::getDx(index_t j) const {
385 213050 : return dx_[j];
386 : }
387 :
388 37 : double GridBase::getBinVolume() const {
389 : double vol=1.;
390 102 : for(unsigned i=0; i<dx_.size(); ++i) {
391 65 : vol*=dx_[i];
392 : }
393 37 : return vol;
394 : }
395 :
396 26440 : std::vector<bool> GridBase::getIsPeriodic() const {
397 26440 : return pbc_;
398 : }
399 :
400 1050474 : std::vector<unsigned> GridBase::getNbin() const {
401 1050474 : return nbin_;
402 : }
403 :
404 9436 : std::vector<std::string> GridBase::getArgNames() const {
405 9436 : return argnames;
406 : }
407 :
408 18384374 : unsigned GridBase::getDimension() const {
409 18384374 : return dimension_;
410 : }
411 :
412 48272 : unsigned GridBase::getSplineNeighbors(const unsigned* indices, std::size_t indices_size, index_t* neighbors, std::size_t neighbors_size)const {
413 48272 : plumed_assert(indices_size==dimension_);
414 48272 : unsigned nneigh=1<<dimension_; // same as unsigned(pow(2.0,dimension_));
415 48272 : plumed_assert(neighbors_size==nneigh);
416 :
417 : unsigned nneighbors = 0;
418 :
419 : std::array<unsigned,maxdim> nindices;
420 : unsigned inind;
421 180804 : for(unsigned int i=0; i<nneigh; ++i) {
422 : unsigned tmp=i;
423 : inind=0;
424 337040 : for(unsigned int j=0; j<dimension_; ++j) {
425 204508 : unsigned i0=tmp%2+indices[j];
426 204508 : tmp/=2;
427 204508 : if(!pbc_[j] && i0==nbin_[j]) {
428 4 : continue;
429 : }
430 204504 : if( pbc_[j] && i0==nbin_[j]) {
431 : i0=0;
432 : }
433 204504 : nindices[inind++]=i0;
434 : }
435 132532 : if(inind==dimension_) {
436 132528 : neighbors[nneighbors++]=getIndex(nindices.data(),dimension_);
437 : }
438 : }
439 48272 : return nneighbors;
440 : }
441 :
442 29207 : std::vector<GridBase::index_t> GridBase::getNearestNeighbors(const index_t index) const {
443 29207 : std::vector<index_t> nearest_neighs = std::vector<index_t>();
444 87620 : for (unsigned i = 0; i < dimension_; i++) {
445 58413 : std::vector<unsigned> neighsneeded = std::vector<unsigned>(dimension_, 0);
446 58413 : neighsneeded[i] = 1;
447 58413 : std::vector<index_t> singledim_nearest_neighs = getNeighbors(index, neighsneeded);
448 232018 : for (unsigned j = 0; j < singledim_nearest_neighs.size(); j++) {
449 173605 : index_t neigh = singledim_nearest_neighs[j];
450 173605 : if (neigh != index) {
451 115192 : nearest_neighs.push_back(neigh);
452 : }
453 : }
454 : }
455 29207 : return nearest_neighs;
456 : }
457 :
458 0 : std::vector<GridBase::index_t> GridBase::getNearestNeighbors(const std::vector<unsigned> &indices) const {
459 : plumed_dbg_assert(indices.size() == dimension_);
460 0 : return getNearestNeighbors(getIndex(indices));
461 : }
462 :
463 0 : void GridBase::addKernel( const KernelFunctions& kernel ) {
464 : plumed_dbg_assert( kernel.ndim()==dimension_ );
465 0 : std::vector<unsigned> nneighb=kernel.getSupport( dx_ );
466 0 : std::vector<index_t> neighbors=getNeighbors( kernel.getCenter(), nneighb );
467 0 : std::vector<double> xx( dimension_ );
468 0 : std::vector<std::unique_ptr<Value>> vv( dimension_ );
469 : std::string str_min, str_max;
470 0 : for(unsigned i=0; i<dimension_; ++i) {
471 0 : vv[i]=Tools::make_unique<Value>();
472 0 : if( pbc_[i] ) {
473 0 : Tools::convert(min_[i],str_min);
474 0 : Tools::convert(max_[i],str_max);
475 0 : vv[i]->setDomain( str_min, str_max );
476 : } else {
477 0 : vv[i]->setNotPeriodic();
478 : }
479 : }
480 :
481 : // vv_ptr contains plain pointers obtained from vv.
482 : // this is the simplest way to replace a unique_ptr here.
483 : // perhaps the interface of kernel.evaluate() should be changed
484 : // in order to accept a std::vector<std::unique_ptr<Value>>
485 0 : auto vv_ptr=Tools::unique2raw(vv);
486 :
487 0 : std::vector<double> der( dimension_ );
488 0 : for(unsigned i=0; i<neighbors.size(); ++i) {
489 0 : index_t ineigh=neighbors[i];
490 0 : getPoint( ineigh, xx );
491 0 : for(unsigned j=0; j<dimension_; ++j) {
492 0 : vv[j]->set(xx[j]);
493 : }
494 0 : double newval = kernel.evaluate( vv_ptr, der, usederiv_ );
495 0 : if( usederiv_ ) {
496 0 : addValueAndDerivatives( ineigh, newval, der );
497 : } else {
498 0 : addValue( ineigh, newval );
499 : }
500 : }
501 0 : }
502 :
503 1193613 : double GridBase::getValue(const std::vector<unsigned> & indices) const {
504 1193613 : return getValue(getIndex(indices));
505 : }
506 :
507 10645 : double GridBase::getValue(const std::vector<double> & x) const {
508 10645 : if(!dospline_) {
509 127 : return getValue(getIndex(x));
510 : } else {
511 10518 : std::vector<double> der(dimension_);
512 10518 : return getValueAndDerivatives(x,der);
513 : }
514 : }
515 :
516 3100944 : double GridBase::getValueAndDerivatives(index_t index, std::vector<double>& der) const {
517 : plumed_dbg_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
518 3100944 : der.resize(dimension_);
519 3100944 : return getValueAndDerivatives(index,der.data(),der.size());
520 : }
521 :
522 2527708 : double GridBase::getValueAndDerivatives(const std::vector<unsigned> & indices, std::vector<double>& der) const {
523 2527708 : return getValueAndDerivatives(getIndex(indices),der);
524 : }
525 :
526 48818 : double GridBase::getValueAndDerivatives(const std::vector<double> & x, std::vector<double>& der) const {
527 : plumed_dbg_assert(der.size()==dimension_ && usederiv_);
528 :
529 48818 : if(dospline_) {
530 : double X,X2,X3,value;
531 : std::array<double,maxdim> fd, C, D;
532 : std::array<double,maxdim> dder;
533 : // reset
534 : value=0.0;
535 114538 : for(unsigned int i=0; i<dimension_; ++i) {
536 66266 : der[i]=0.0;
537 : }
538 :
539 : std::array<unsigned,maxdim> indices;
540 48272 : getIndices(x, indices.data(),dimension_);
541 : std::array<double,maxdim> xfloor;
542 48272 : getPoint(indices.data(), dimension_, xfloor.data(),dimension_);
543 48272 : gch::small_vector<index_t,16> neigh(1<<dimension_); // pow(2,dimension_); up to dimension 4 will stay on stack
544 48272 : auto nneigh = getSplineNeighbors(indices.data(),dimension_, neigh.data(), neigh.size());
545 :
546 : // loop over neighbors
547 : std::array<unsigned,maxdim> nindices;
548 180800 : for(unsigned int ipoint=0; ipoint<nneigh; ++ipoint) {
549 132528 : double grid=getValueAndDerivatives(neigh[ipoint],dder.data(),dimension_);
550 132528 : getIndices(neigh[ipoint], nindices.data(), dimension_);
551 : double ff=1.0;
552 :
553 337032 : for(unsigned j=0; j<dimension_; ++j) {
554 : int x0=1;
555 204504 : if(nindices[j]==indices[j]) {
556 : x0=0;
557 : }
558 204504 : double dx=getDx(j);
559 204504 : X=std::abs((x[j]-xfloor[j])/dx-(double)x0);
560 204504 : X2=X*X;
561 204504 : X3=X2*X;
562 : double yy;
563 204504 : if(std::abs(grid)<0.0000001) {
564 : yy=0.0;
565 : } else {
566 200190 : yy=-dder[j]/grid;
567 : }
568 204504 : C[j]=(1.0-3.0*X2+2.0*X3) - (x0?-1.0:1.0)*yy*(X-2.0*X2+X3)*dx;
569 204504 : D[j]=( -6.0*X +6.0*X2) - (x0?-1.0:1.0)*yy*(1.0-4.0*X +3.0*X2)*dx;
570 204504 : D[j]*=(x0?-1.0:1.0)/dx;
571 204504 : ff*=C[j];
572 : }
573 337032 : for(unsigned j=0; j<dimension_; ++j) {
574 204504 : fd[j]=D[j];
575 552960 : for(unsigned i=0; i<dimension_; ++i)
576 348456 : if(i!=j) {
577 143952 : fd[j]*=C[i];
578 : }
579 : }
580 132528 : value+=grid*ff;
581 337032 : for(unsigned j=0; j<dimension_; ++j) {
582 204504 : der[j]+=grid*fd[j];
583 : }
584 : }
585 : return value;
586 : } else {
587 546 : return getValueAndDerivatives(getIndex(x),der);
588 : }
589 : }
590 :
591 0 : void GridBase::setValue(const std::vector<unsigned> & indices, double value) {
592 0 : setValue(getIndex(indices),value);
593 0 : }
594 :
595 0 : void GridBase::setValueAndDerivatives(const std::vector<unsigned> & indices, double value, std::vector<double>& der) {
596 0 : setValueAndDerivatives(getIndex(indices),value,der);
597 0 : }
598 :
599 0 : void GridBase::addValue(const std::vector<unsigned> & indices, double value) {
600 0 : addValue(getIndex(indices),value);
601 0 : }
602 :
603 0 : void GridBase::addValueAndDerivatives(const std::vector<unsigned> & indices, double value, std::vector<double>& der) {
604 0 : addValueAndDerivatives(getIndex(indices),value,der);
605 0 : }
606 :
607 0 : void GridBase::addValueAndDerivatives(index_t index, double value, const double* der, std::size_t der_size) {
608 0 : std::vector<double> d(der,der+der_size);
609 0 : addValueAndDerivatives(index,value,d);
610 0 : }
611 :
612 1180 : void GridBase::writeHeader(OFile& ofile) {
613 2685 : for(unsigned i=0; i<dimension_; ++i) {
614 3010 : ofile.addConstantField("min_" + argnames[i]);
615 3010 : ofile.addConstantField("max_" + argnames[i]);
616 3010 : ofile.addConstantField("nbins_" + argnames[i]);
617 3010 : ofile.addConstantField("periodic_" + argnames[i]);
618 : }
619 1180 : }
620 :
621 1 : void Grid::clear() {
622 1 : grid_.assign(maxsize_,0.0);
623 1 : if(usederiv_) {
624 0 : der_.assign(maxsize_*dimension_,0.0);
625 : }
626 1 : }
627 :
628 1180 : void Grid::writeToFile(OFile& ofile) {
629 1180 : std::vector<double> xx(dimension_);
630 1180 : std::vector<double> der(dimension_);
631 : double f;
632 1180 : writeHeader(ofile);
633 2562753 : for(index_t i=0; i<getSize(); ++i) {
634 2561573 : xx=getPoint(i);
635 2561573 : if(usederiv_) {
636 571866 : f=getValueAndDerivatives(i,der);
637 : } else {
638 1989707 : f=getValue(i);
639 : }
640 4971197 : if(i>0 && dimension_>1 && getIndices(i)[dimension_-2]==0) {
641 25000 : ofile.printf("\n");
642 : }
643 7570135 : for(unsigned j=0; j<dimension_; ++j) {
644 10017124 : ofile.printField("min_" + argnames[j], str_min_[j] );
645 10017124 : ofile.printField("max_" + argnames[j], str_max_[j] );
646 10017124 : ofile.printField("nbins_" + argnames[j], static_cast<int>(nbin_[j]) );
647 5008562 : if( pbc_[j] ) {
648 4101894 : ofile.printField("periodic_" + argnames[j], "true" );
649 : } else {
650 5915230 : ofile.printField("periodic_" + argnames[j], "false" );
651 : }
652 : }
653 7570135 : for(unsigned j=0; j<dimension_; ++j) {
654 5008562 : ofile.fmtField(" "+fmt_);
655 5008562 : ofile.printField(argnames[j],xx[j]);
656 : }
657 2561573 : ofile.fmtField(" "+fmt_);
658 2561573 : ofile.printField(funcname,f);
659 2561573 : if(usederiv_)
660 1687906 : for(unsigned j=0; j<dimension_; ++j) {
661 1116040 : ofile.fmtField(" "+fmt_);
662 2232080 : ofile.printField("der_" + argnames[j],der[j]);
663 : }
664 2561573 : ofile.printField();
665 : }
666 1180 : }
667 :
668 0 : void GridBase::writeCubeFile(OFile& ofile, const double& lunit) {
669 0 : plumed_assert( dimension_==3 );
670 0 : ofile.printf("PLUMED CUBE FILE\n");
671 0 : ofile.printf("OUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z\n");
672 : // Number of atoms followed by position of origin (origin set so that center of grid is in center of cell)
673 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]));
674 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
675 0 : ofile.printf("%u %f %f %f\n",nbin_[1],0.0,lunit*dx_[1],0.0); // shape of voxel
676 0 : ofile.printf("%u %f %f %f\n",nbin_[2],0.0,0.0,lunit*dx_[2]);
677 0 : ofile.printf("%d %f %f %f\n",1,0.0,0.0,0.0); // Fake atom otherwise VMD doesn't work
678 0 : std::vector<unsigned> pp(3);
679 0 : for(pp[0]=0; pp[0]<nbin_[0]; ++pp[0]) {
680 0 : for(pp[1]=0; pp[1]<nbin_[1]; ++pp[1]) {
681 0 : for(pp[2]=0; pp[2]<nbin_[2]; ++pp[2]) {
682 0 : ofile.printf("%f ",getValue(pp) );
683 0 : if(pp[2]%6==5) {
684 0 : ofile.printf("\n");
685 : }
686 : }
687 0 : ofile.printf("\n");
688 : }
689 : }
690 0 : }
691 :
692 0 : std::unique_ptr<GridBase> GridBase::create(const std::string& funcl, const std::vector<Value*> & args, IFile& ifile,
693 : const std::vector<std::string> & gmin,const std::vector<std::string> & gmax,
694 : const std::vector<unsigned> & nbin,bool dosparse, bool dospline, bool doder) {
695 0 : std::unique_ptr<GridBase> grid=GridBase::create(funcl,args,ifile,dosparse,dospline,doder);
696 0 : std::vector<unsigned> cbin( grid->getNbin() );
697 0 : std::vector<std::string> cmin( grid->getMin() ), cmax( grid->getMax() );
698 0 : for(unsigned i=0; i<args.size(); ++i) {
699 0 : plumed_massert( cmin[i]==gmin[i], "mismatched grid min" );
700 0 : plumed_massert( cmax[i]==gmax[i], "mismatched grid max" );
701 0 : if( args[i]->isPeriodic() ) {
702 0 : plumed_massert( cbin[i]==nbin[i], "mismatched grid nbins" );
703 : } else {
704 0 : plumed_massert( (cbin[i]-1)==nbin[i], "mismatched grid nbins");
705 : }
706 : }
707 0 : return grid;
708 0 : }
709 :
710 93 : std::unique_ptr<GridBase> GridBase::create(const std::string& funcl, const std::vector<Value*> & args, IFile& ifile, bool dosparse, bool dospline, bool doder) {
711 93 : std::unique_ptr<GridBase> grid;
712 93 : unsigned nvar=args.size();
713 : bool hasder=false;
714 : std::string pstring;
715 94 : std::vector<int> gbin1(nvar);
716 94 : std::vector<unsigned> gbin(nvar);
717 94 : std::vector<std::string> labels(nvar),gmin(nvar),gmax(nvar);
718 93 : std::vector<std::string> fieldnames;
719 93 : ifile.scanFieldList( fieldnames );
720 :
721 : // Retrieve names for fields
722 214 : for(unsigned i=0; i<args.size(); ++i) {
723 121 : labels[i]=args[i]->getName();
724 : bool found=false;
725 149 : for(unsigned j=0; j<fieldnames.size(); ++j) {
726 149 : if( labels[i]==fieldnames[j] ) {
727 : found=true;
728 : break;
729 : }
730 : }
731 : // This is a change to ensure that we can deal with old style names for multicolvars
732 121 : std::size_t und = labels[i].find_first_of("_");
733 121 : if( !found && und!=std::string::npos ) {
734 0 : labels[i] = labels[i].substr(0,und) + "." + labels[i].substr(und+1);
735 0 : for(unsigned j=0; j<fieldnames.size(); ++j) {
736 0 : if( labels[i]==fieldnames[j] ) {
737 : found=true;
738 : break;
739 : }
740 : }
741 : }
742 : }
743 : // And read the stuff from the header
744 93 : if( !ifile.FieldExist( funcl ) ) {
745 : std::string fields;
746 8 : for(const auto & f : fieldnames) {
747 15 : fields += " " + f;
748 : }
749 5 : plumed_merror("there is no column labelled " + funcl + " in the grid file " + ifile.getPath() + " (its columns are:" + fields + "). "
750 : "The column is named after the label of the action that wrote the grid: check that the action has the same label");
751 : }
752 212 : for(unsigned i=0; i<args.size(); ++i) {
753 240 : ifile.scanField( "min_" + labels[i], gmin[i]);
754 240 : ifile.scanField( "max_" + labels[i], gmax[i]);
755 240 : ifile.scanField( "periodic_" + labels[i], pstring );
756 240 : ifile.scanField( "nbins_" + labels[i], gbin1[i]);
757 120 : plumed_assert( gbin1[i]>0 );
758 120 : if( args[i]->isPeriodic() ) {
759 45 : plumed_massert( pstring=="true", "input value is periodic but grid is not");
760 : std::string pmin, pmax;
761 45 : args[i]->getDomain( pmin, pmax );
762 45 : gbin[i]=gbin1[i];
763 45 : if( pmin!=gmin[i] || pmax!=gmax[i] ) {
764 0 : plumed_merror("mismatch between grid boundaries and periods of values");
765 : }
766 : } else {
767 75 : gbin[i]=gbin1[i]-1; // Note header in grid file indicates one more bin that there should be when data is not periodic
768 75 : plumed_massert( pstring=="false", "input value is not periodic but grid is");
769 : }
770 120 : hasder=ifile.FieldExist( "der_" + labels[i] );
771 120 : if( doder && !hasder ) {
772 0 : plumed_merror("missing derivatives from grid file");
773 : }
774 148 : for(unsigned j=0; j<fieldnames.size(); ++j) {
775 176 : for(unsigned k=i+1; k<args.size(); ++k) {
776 28 : if( fieldnames[j]==labels[k] ) {
777 0 : plumed_merror("arguments in input are not in same order as in grid file");
778 : }
779 : }
780 148 : if( fieldnames[j]==labels[i] ) {
781 : break;
782 : }
783 : }
784 : }
785 :
786 92 : if(!dosparse) {
787 91 : grid=Tools::make_unique<Grid>(funcl,args,gmin,gmax,gbin,dospline,doder);
788 : } else {
789 1 : grid=Tools::make_unique<SparseGrid>(funcl,args,gmin,gmax,gbin,dospline,doder);
790 : }
791 :
792 92 : std::vector<double> xx(nvar),dder(nvar);
793 92 : std::vector<double> dx=grid->getDx();
794 : double f,x;
795 1138531 : while( ifile.scanField(funcl,f) ) {
796 3410248 : for(unsigned i=0; i<nvar; ++i) {
797 2271809 : ifile.scanField(labels[i],x);
798 2271809 : xx[i]=x+dx[i]/2.0;
799 4543618 : ifile.scanField( "min_" + labels[i], gmin[i]);
800 4543618 : ifile.scanField( "max_" + labels[i], gmax[i]);
801 4543618 : ifile.scanField( "nbins_" + labels[i], gbin1[i]);
802 4543618 : ifile.scanField( "periodic_" + labels[i], pstring );
803 : }
804 1138439 : if(hasder) {
805 3284226 : for(unsigned i=0; i<nvar; ++i) {
806 4376792 : ifile.scanField( "der_" + labels[i], dder[i] );
807 : }
808 : }
809 1138439 : index_t index=grid->getIndex(xx);
810 1138439 : if(doder) {
811 1095722 : grid->setValueAndDerivatives(index,f,dder);
812 : } else {
813 42717 : grid->setValue(index,f);
814 : }
815 1138439 : ifile.scanField();
816 : }
817 92 : return grid;
818 97 : }
819 :
820 861 : double Grid::getMinValue() const {
821 : double minval;
822 : minval=DBL_MAX;
823 2258311 : for(index_t i=0; i<grid_.size(); ++i) {
824 2257450 : if(grid_[i]<minval) {
825 : minval=grid_[i];
826 : }
827 : }
828 861 : return minval;
829 : }
830 :
831 1202 : double Grid::getMaxValue() const {
832 : double maxval;
833 : maxval=DBL_MIN;
834 11595330 : for(index_t i=0; i<grid_.size(); ++i) {
835 11594128 : if(grid_[i]>maxval) {
836 : maxval=grid_[i];
837 : }
838 : }
839 1202 : return maxval;
840 : }
841 :
842 607 : void Grid::scaleAllValuesAndDerivatives( const double& scalef ) {
843 607 : if(usederiv_) {
844 33548 : for(index_t i=0; i<grid_.size(); ++i) {
845 33524 : grid_[i]*=scalef;
846 99570 : for(unsigned j=0; j<dimension_; ++j) {
847 66046 : der_[i*dimension_+j]*=scalef;
848 : }
849 : }
850 : } else {
851 1572834 : for(index_t i=0; i<grid_.size(); ++i) {
852 1572251 : grid_[i]*=scalef;
853 : }
854 : }
855 607 : }
856 :
857 0 : void Grid::logAllValuesAndDerivatives( const double& scalef ) {
858 0 : if(usederiv_) {
859 0 : for(index_t i=0; i<grid_.size(); ++i) {
860 0 : grid_[i] = scalef*std::log(grid_[i]);
861 0 : for(unsigned j=0; j<dimension_; ++j) {
862 0 : der_[i*dimension_+j] = scalef/der_[i*dimension_+j];
863 : }
864 : }
865 : } else {
866 0 : for(index_t i=0; i<grid_.size(); ++i) {
867 0 : grid_[i] = scalef*std::log(grid_[i]);
868 : }
869 : }
870 0 : }
871 :
872 1329 : void Grid::setMinToZero() {
873 1329 : double min=grid_[0];
874 3518541 : for(index_t i=1; i<grid_.size(); ++i)
875 3517212 : if(grid_[i]<min) {
876 : min=grid_[i];
877 : }
878 3519870 : for(index_t i=0; i<grid_.size(); ++i) {
879 3518541 : grid_[i] -= min;
880 : }
881 1329 : }
882 :
883 1 : void Grid::applyFunctionAllValuesAndDerivatives( double (*func)(double val), double (*funcder)(double valder) ) {
884 1 : if(usederiv_) {
885 901 : for(index_t i=0; i<grid_.size(); ++i) {
886 900 : grid_[i]=func(grid_[i]);
887 2700 : for(unsigned j=0; j<dimension_; ++j) {
888 1800 : der_[i*dimension_+j]=funcder(der_[i*dimension_+j]);
889 : }
890 : }
891 : } else {
892 0 : for(index_t i=0; i<grid_.size(); ++i) {
893 0 : grid_[i]=func(grid_[i]);
894 : }
895 : }
896 1 : }
897 :
898 0 : double Grid::getDifferenceFromContour( const std::vector<double>& x, std::vector<double>& der ) const {
899 0 : return getValueAndDerivatives( x, der ) - contour_location;
900 : }
901 :
902 0 : void Grid::findSetOfPointsOnContour(const double& target, const std::vector<bool>& nosearch,
903 : unsigned& npoints, std::vector<std::vector<double> >& points ) {
904 : // Set contour location for function
905 0 : contour_location=target;
906 : // Resize points to maximum possible value
907 0 : points.resize( dimension_*maxsize_ );
908 :
909 : // Two points for search
910 0 : std::vector<unsigned> ind(dimension_);
911 0 : std::vector<double> direction( dimension_, 0 );
912 :
913 : // Run over whole grid
914 0 : npoints=0;
915 : RootFindingBase<Grid> mymin( this );
916 0 : for(unsigned i=0; i<maxsize_; ++i) {
917 0 : for(unsigned j=0; j<dimension_; ++j) {
918 0 : ind[j]=getIndices(i)[j];
919 : }
920 :
921 : // Get the value of a point on the grid
922 0 : double val1=getValue(i) - target;
923 :
924 : // Now search for contour in each direction
925 : bool edge=false;
926 0 : for(unsigned j=0; j<dimension_; ++j) {
927 0 : if( nosearch[j] ) {
928 0 : continue ;
929 : }
930 : // Make sure we don't search at the edge of the grid
931 0 : if( !pbc_[j] && (ind[j]+1)==nbin_[j] ) {
932 0 : continue;
933 0 : } else if( (ind[j]+1)==nbin_[j] ) {
934 : edge=true;
935 0 : ind[j]=0;
936 : } else {
937 0 : ind[j]+=1;
938 : }
939 0 : double val2=getValue(ind) - target;
940 0 : if( val1*val2<0 ) {
941 : // Use initial point location as first guess for search
942 0 : points[npoints].resize(dimension_);
943 0 : for(unsigned k=0; k<dimension_; ++k) {
944 0 : points[npoints][k]=getPoint(i)[k];
945 : }
946 : // Setup direction vector
947 0 : direction[j]=0.999999999*dx_[j];
948 : // And do proper search for contour point
949 0 : mymin.linesearch( direction, points[npoints], &Grid::getDifferenceFromContour );
950 0 : direction[j]=0.0;
951 0 : npoints++;
952 : }
953 0 : if( pbc_[j] && edge ) {
954 : edge=false;
955 0 : ind[j]=nbin_[j]-1;
956 : } else {
957 0 : ind[j]-=1;
958 : }
959 : }
960 : }
961 0 : }
962 :
963 : /// OVERRIDES ARE BELOW
964 :
965 27984931 : Grid::index_t Grid::getSize() const {
966 27984931 : return maxsize_;
967 : }
968 :
969 35819167 : double Grid::getValue(index_t index) const {
970 : plumed_dbg_assert(index<maxsize_);
971 35819167 : return grid_[index];
972 : }
973 :
974 3215376 : double Grid::getValueAndDerivatives(index_t index, double* der,std::size_t der_size) const {
975 : plumed_dbg_assert(index<maxsize_ && usederiv_ && der_size==dimension_);
976 7034252 : for(unsigned i=0; i<dimension_; i++) {
977 3818876 : der[i]=der_[dimension_*index+i];
978 : }
979 3215376 : return grid_[index];
980 : }
981 :
982 6444104 : void Grid::setValue(index_t index, double value) {
983 : plumed_dbg_assert(index<maxsize_ && !usederiv_);
984 6444104 : grid_[index]=value;
985 6444104 : }
986 :
987 2592840 : void Grid::setValueAndDerivatives(index_t index, double value, std::vector<double>& der) {
988 : plumed_dbg_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
989 2592840 : grid_[index]=value;
990 7727932 : for(unsigned i=0; i<dimension_; i++) {
991 5135092 : der_[dimension_*index+i]=der[i];
992 : }
993 2592840 : }
994 :
995 1 : void Grid::addValue(index_t index, double value) {
996 : plumed_dbg_assert(index<maxsize_ && !usederiv_);
997 1 : grid_[index]+=value;
998 1 : }
999 :
1000 1058017 : void Grid::addValueAndDerivatives(index_t index, double value, std::vector<double>& der) {
1001 : plumed_dbg_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
1002 1058017 : grid_[index]+=value;
1003 3086538 : for(unsigned int i=0; i<dimension_; ++i) {
1004 2028521 : der_[index*dimension_+i]+=der[i];
1005 : }
1006 1058017 : }
1007 :
1008 4453920 : void Grid::addValueAndDerivatives(index_t index, double value, const double* der, std::size_t der_size) {
1009 : plumed_dbg_assert(index<maxsize_ && usederiv_ && der_size==dimension_);
1010 4453920 : grid_[index]+=value;
1011 13095620 : for(unsigned int i=0; i<dimension_; ++i) {
1012 8641700 : der_[index*dimension_+i]+=der[i];
1013 : }
1014 4453920 : }
1015 :
1016 123 : void Grid::resizeData(const std::vector<unsigned> & oldnbin, const std::vector<long long> & shift) {
1017 123 : std::vector<double> newgrid(maxsize_,0.0);
1018 123 : std::vector<double> newder;
1019 123 : if(usederiv_) {
1020 123 : newder.assign(maxsize_*dimension_,0.0);
1021 : }
1022 : index_t oldsize=1;
1023 291 : for(unsigned i=0; i<dimension_; ++i) {
1024 168 : oldsize*=oldnbin[i];
1025 : }
1026 : // loop over the old points in column-major order
1027 123 : std::vector<unsigned> k(dimension_,0);
1028 152290 : for(index_t old=0; old<oldsize; ++old) {
1029 : bool kept=true;
1030 : index_t index=0;
1031 446599 : for(unsigned i=dimension_; i-->0;) {
1032 295714 : const long long knew=static_cast<long long>(k[i])+shift[i];
1033 295714 : if(knew<0 || knew>=static_cast<long long>(nbin_[i])) {
1034 : kept=false;
1035 : break;
1036 : }
1037 294432 : index=index*nbin_[i]+static_cast<index_t>(knew);
1038 : }
1039 152167 : if(kept) {
1040 150885 : newgrid[index]=grid_[old];
1041 150885 : if(usederiv_) {
1042 445317 : for(unsigned i=0; i<dimension_; ++i) {
1043 294432 : newder[index*dimension_+i]=der_[old*dimension_+i];
1044 : }
1045 : }
1046 : }
1047 153572 : for(unsigned i=0; i<dimension_; ++i) {
1048 153449 : if(++k[i]<oldnbin[i]) {
1049 : break;
1050 : }
1051 1405 : k[i]=0;
1052 : }
1053 : }
1054 : grid_.swap(newgrid);
1055 : der_.swap(newder);
1056 123 : }
1057 :
1058 9 : void SparseGrid::resizeData(const std::vector<unsigned> & oldnbin, const std::vector<long long> & shift) {
1059 : // new index of the point with index old before resizing, false if it was removed
1060 17190 : auto newIndex=[&](index_t old, index_t & index) {
1061 17190 : std::vector<long long> k(dimension_);
1062 51570 : for(unsigned i=0; i<dimension_; ++i) {
1063 34380 : k[i]=static_cast<long long>(old%oldnbin[i])+shift[i];
1064 34380 : old/=oldnbin[i];
1065 34380 : if(k[i]<0 || k[i]>=static_cast<long long>(nbin_[i])) {
1066 : return false;
1067 : }
1068 : }
1069 17190 : index=0;
1070 51570 : for(unsigned i=dimension_; i-->0;) {
1071 34380 : index=index*nbin_[i]+static_cast<index_t>(k[i]);
1072 : }
1073 : return true;
1074 9 : };
1075 : std::map<index_t,double> newmap;
1076 8604 : for(const auto & it : map_) {
1077 : index_t index;
1078 8595 : if(newIndex(it.first,index)) {
1079 8595 : newmap[index]=it.second;
1080 : }
1081 : }
1082 : std::map<index_t,std::vector<double>> newder;
1083 8604 : for(auto & it : der_) {
1084 : index_t index;
1085 8595 : if(newIndex(it.first,index)) {
1086 8595 : newder[index].swap(it.second);
1087 : }
1088 : }
1089 : map_.swap(newmap);
1090 : der_.swap(newder);
1091 9 : }
1092 :
1093 0 : Grid::index_t SparseGrid::getSize() const {
1094 0 : return map_.size();
1095 : }
1096 :
1097 0 : Grid::index_t SparseGrid::getMaxSize() const {
1098 0 : return maxsize_;
1099 : }
1100 :
1101 0 : double SparseGrid::getValue(index_t index)const {
1102 0 : plumed_assert(index<maxsize_);
1103 : double value=0.0;
1104 : const auto it=map_.find(index);
1105 0 : if(it!=map_.end()) {
1106 0 : value=it->second;
1107 : }
1108 0 : return value;
1109 : }
1110 :
1111 18456 : double SparseGrid::getValueAndDerivatives(index_t index, double* der, std::size_t der_size)const {
1112 18456 : plumed_assert(index<maxsize_ && usederiv_ && der_size==dimension_);
1113 : double value=0.0;
1114 50108 : for(unsigned int i=0; i<dimension_; ++i) {
1115 31652 : der[i]=0.0;
1116 : }
1117 : const auto it=map_.find(index);
1118 18456 : if(it!=map_.end()) {
1119 18290 : value=it->second;
1120 : }
1121 : const auto itder=der_.find(index);
1122 18456 : if(itder!=der_.end()) {
1123 : const auto & second(itder->second);
1124 49660 : for(unsigned i=0; i<second.size(); i++) {
1125 31370 : der[i]=itder->second[i];
1126 : }
1127 : }
1128 18456 : return value;
1129 : }
1130 :
1131 0 : void SparseGrid::setValue(index_t index, double value) {
1132 0 : plumed_assert(index<maxsize_ && !usederiv_);
1133 0 : map_[index]=value;
1134 0 : }
1135 :
1136 37 : void SparseGrid::setValueAndDerivatives(index_t index, double value, std::vector<double>& der) {
1137 37 : plumed_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
1138 37 : map_[index]=value;
1139 37 : der_[index]=der;
1140 37 : }
1141 :
1142 0 : void SparseGrid::addValue(index_t index, double value) {
1143 0 : plumed_assert(index<maxsize_ && !usederiv_);
1144 0 : map_[index]+=value;
1145 0 : }
1146 :
1147 0 : void SparseGrid::addValueAndDerivatives(index_t index, double value, std::vector<double>& der) {
1148 0 : plumed_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
1149 0 : map_[index]+=value;
1150 0 : der_[index].resize(dimension_);
1151 0 : for(unsigned int i=0; i<dimension_; ++i) {
1152 0 : der_[index][i]+=der[i];
1153 : }
1154 0 : }
1155 :
1156 796842 : void SparseGrid::addValueAndDerivatives(index_t index, double value, const double* der, std::size_t der_size) {
1157 796842 : plumed_assert(index<maxsize_ && usederiv_ && der_size==dimension_);
1158 796842 : map_[index]+=value;
1159 796842 : der_[index].resize(dimension_);
1160 2373643 : for(unsigned int i=0; i<dimension_; ++i) {
1161 1576801 : der_[index][i]+=der[i];
1162 : }
1163 796842 : }
1164 :
1165 0 : void SparseGrid::writeToFile(OFile& ofile) {
1166 0 : std::vector<double> xx(dimension_);
1167 0 : std::vector<double> der(dimension_);
1168 : double f;
1169 0 : writeHeader(ofile);
1170 0 : ofile.fmtField(" "+fmt_);
1171 0 : for(const auto & it : map_) {
1172 0 : index_t i=it.first;
1173 0 : xx=getPoint(i);
1174 0 : if(usederiv_) {
1175 0 : f=getValueAndDerivatives(i,der);
1176 : } else {
1177 0 : f=getValue(i);
1178 : }
1179 0 : if(i>0 && dimension_>1 && getIndices(i)[dimension_-2]==0) {
1180 0 : ofile.printf("\n");
1181 : }
1182 0 : for(unsigned j=0; j<dimension_; ++j) {
1183 0 : ofile.printField("min_" + argnames[j], str_min_[j] );
1184 0 : ofile.printField("max_" + argnames[j], str_max_[j] );
1185 0 : ofile.printField("nbins_" + argnames[j], static_cast<int>(nbin_[j]) );
1186 0 : if( pbc_[j] ) {
1187 0 : ofile.printField("periodic_" + argnames[j], "true" );
1188 : } else {
1189 0 : ofile.printField("periodic_" + argnames[j], "false" );
1190 : }
1191 : }
1192 0 : for(unsigned j=0; j<dimension_; ++j) {
1193 0 : ofile.printField(argnames[j],xx[j]);
1194 : }
1195 0 : ofile.printField(funcname, f);
1196 0 : if(usederiv_) {
1197 0 : for(unsigned j=0; j<dimension_; ++j) {
1198 0 : ofile.printField("der_" + argnames[j],der[j]);
1199 : }
1200 : }
1201 0 : ofile.printField();
1202 : }
1203 0 : }
1204 :
1205 0 : double SparseGrid::getMinValue() const {
1206 : double minval;
1207 : minval=0.0;
1208 0 : for(auto const & i : map_) {
1209 0 : if(i.second<minval) {
1210 : minval=i.second;
1211 : }
1212 : }
1213 0 : return minval;
1214 : }
1215 :
1216 445 : double SparseGrid::getMaxValue() const {
1217 : double maxval;
1218 : maxval=0.0;
1219 1095898 : for(auto const & i : map_) {
1220 1095453 : if(i.second>maxval) {
1221 : maxval=i.second;
1222 : }
1223 : }
1224 445 : return maxval;
1225 : }
1226 :
1227 1010062 : void Grid::projectOnLowDimension(double &val, std::vector<int> &vHigh, WeightBase * ptr2obj ) {
1228 : unsigned i=0;
1229 3016945 : for(i=0; i<vHigh.size(); i++) {
1230 2015726 : if(vHigh[i]<0) { // this bin needs to be integrated out
1231 : // parallelize here???
1232 1010062 : for(unsigned j=0; j<(getNbin())[i]; j++) {
1233 1001219 : vHigh[i]=int(j);
1234 1001219 : projectOnLowDimension(val,vHigh,ptr2obj); // recursive function: this is the core of the mechanism
1235 1001219 : vHigh[i]=-1;
1236 : }
1237 : return; //
1238 : }
1239 : }
1240 : // when there are no more bin to dig in then retrieve the value
1241 1001219 : if(i==vHigh.size()) {
1242 1001219 : std::vector<unsigned> vv(vHigh.size());
1243 3003657 : for(unsigned j=0; j<vHigh.size(); j++) {
1244 2002438 : vv[j]=unsigned(vHigh[j]);
1245 : }
1246 : //
1247 : // this is the real assignment !!!!! (hack this to have bias or other stuff)
1248 : //
1249 :
1250 : // this case: produce fes
1251 : //val+=exp(beta*getValue(vv)) ;
1252 1001219 : double myv=getValue(vv);
1253 1001219 : val=ptr2obj->projectInnerLoop(val,myv) ;
1254 : // to be added: bias (same as before without negative sign)
1255 : //std::cerr<<" VAL: "<<val <<endl;
1256 : }
1257 : }
1258 :
1259 81 : Grid Grid::project(const std::vector<std::string> & proj, WeightBase *ptr2obj ) {
1260 : // find extrema only for the projection
1261 81 : std::vector<std::string> smallMin,smallMax;
1262 81 : std::vector<unsigned> smallBin;
1263 81 : std::vector<unsigned> dimMapping;
1264 81 : std::vector<bool> smallIsPeriodic;
1265 81 : std::vector<std::string> smallName;
1266 :
1267 : // check if the two key methods are there
1268 : WeightBase* pp = dynamic_cast<WeightBase*>(ptr2obj);
1269 81 : if (!pp) {
1270 0 : plumed_merror("This WeightBase is not complete: you need a projectInnerLoop and projectOuterLoop ");
1271 : }
1272 :
1273 162 : for(unsigned j=0; j<proj.size(); j++) {
1274 120 : for(unsigned i=0; i<getArgNames().size(); i++) {
1275 120 : if(proj[j]==getArgNames()[i]) {
1276 : unsigned offset;
1277 : // note that at sizetime the non periodic dimension get a bin more
1278 162 : if(getIsPeriodic()[i]) {
1279 : offset=0;
1280 : } else {
1281 : offset=1;
1282 : }
1283 81 : smallMax.push_back(getMax()[i]);
1284 81 : smallMin.push_back(getMin()[i]);
1285 81 : smallBin.push_back(getNbin()[i]-offset);
1286 81 : smallIsPeriodic.push_back(getIsPeriodic()[i]);
1287 81 : dimMapping.push_back(i);
1288 81 : smallName.push_back(getArgNames()[i]);
1289 81 : break;
1290 : }
1291 : }
1292 : }
1293 81 : Grid smallgrid("projection",smallName,smallMin,smallMax,smallBin,false,false,smallIsPeriodic,smallMin,smallMax);
1294 : // check that the two grids are commensurate
1295 162 : for(unsigned i=0; i<dimMapping.size(); i++) {
1296 81 : plumed_massert( (smallgrid.getMax())[i] == (getMax())[dimMapping[i]], "the two input grids are not compatible in max" );
1297 81 : plumed_massert( (smallgrid.getMin())[i] == (getMin())[dimMapping[i]], "the two input grids are not compatible in min" );
1298 81 : plumed_massert( (smallgrid.getNbin())[i]== (getNbin())[dimMapping[i]], "the two input grids are not compatible in bin" );
1299 : }
1300 81 : std::vector<unsigned> toBeIntegrated;
1301 243 : for(unsigned i=0; i<getArgNames().size(); i++) {
1302 : bool doappend=true;
1303 243 : for(unsigned j=0; j<dimMapping.size(); j++) {
1304 162 : if(dimMapping[j]==i) {
1305 : doappend=false;
1306 : break;
1307 : }
1308 : }
1309 162 : if(doappend) {
1310 81 : toBeIntegrated.push_back(i);
1311 : }
1312 : }
1313 :
1314 : // loop over all the points in the Grid, find the corresponding fixed index, rotate over all the other ones
1315 8924 : for(unsigned i=0; i<smallgrid.getSize(); i++) {
1316 8843 : std::vector<unsigned> v;
1317 8843 : v=smallgrid.getIndices(i);
1318 8843 : std::vector<int> vHigh((getArgNames()).size(),-1);
1319 17686 : for(unsigned j=0; j<dimMapping.size(); j++) {
1320 8843 : vHigh[dimMapping[j]]=int(v[j]);
1321 : }
1322 : // the vector vhigh now contains at the beginning the index of the low dimension and -1 for the values that need to be integrated
1323 8843 : double initval=0.;
1324 8843 : projectOnLowDimension(initval,vHigh, ptr2obj);
1325 8843 : smallgrid.setValue(i,ptr2obj->projectOuterLoop(initval));
1326 : }
1327 :
1328 81 : return smallgrid;
1329 162 : }
1330 :
1331 0 : double Grid::integrate( std::vector<unsigned>& npoints ) {
1332 : plumed_dbg_assert( npoints.size()==dimension_ );
1333 0 : plumed_assert( dospline_ );
1334 :
1335 : unsigned ntotgrid=1;
1336 : double box_vol=1.0;
1337 0 : std::vector<double> ispacing( npoints.size() );
1338 0 : for(unsigned j=0; j<dimension_; ++j) {
1339 0 : if( !pbc_[j] ) {
1340 0 : ispacing[j] = ( max_[j] - dx_[j] - min_[j] ) / static_cast<double>( npoints[j] );
1341 0 : npoints[j]+=1;
1342 : } else {
1343 0 : ispacing[j] = ( max_[j] - min_[j] ) / static_cast<double>( npoints[j] );
1344 : }
1345 0 : ntotgrid*=npoints[j];
1346 0 : box_vol*=ispacing[j];
1347 : }
1348 :
1349 0 : std::vector<double> vals( dimension_ );
1350 0 : std::vector<unsigned> t_index( dimension_ );
1351 : double integral=0.0;
1352 0 : for(unsigned i=0; i<ntotgrid; ++i) {
1353 0 : t_index[0]=(i%npoints[0]);
1354 : unsigned kk=i;
1355 0 : for(unsigned j=1; j<dimension_-1; ++j) {
1356 0 : kk=(kk-t_index[j-1])/npoints[j-1];
1357 0 : t_index[j]=(kk%npoints[j]);
1358 : }
1359 0 : if( dimension_>=2 ) {
1360 0 : t_index[dimension_-1]=((kk-t_index[dimension_-2])/npoints[dimension_-2]);
1361 : }
1362 :
1363 0 : for(unsigned j=0; j<dimension_; ++j) {
1364 0 : vals[j]=min_[j] + t_index[j]*ispacing[j];
1365 : }
1366 :
1367 0 : integral += getValue( vals );
1368 : }
1369 :
1370 0 : return box_vol*integral;
1371 : }
1372 :
1373 0 : void Grid::mpiSumValuesAndDerivatives( Communicator& comm ) {
1374 0 : comm.Sum( grid_ );
1375 0 : for(unsigned i=0; i<der_.size(); ++i) {
1376 0 : comm.Sum( der_[i] );
1377 : }
1378 0 : }
1379 :
1380 :
1381 213811 : bool indexed_lt(std::pair<Grid::index_t, double> const &x, std::pair<Grid::index_t, double> const &y) {
1382 213811 : return x.second < y.second;
1383 : }
1384 :
1385 491 : double GridBase::findMaximalPathMinimum(const std::vector<double> &source, const std::vector<double> &sink) {
1386 : plumed_dbg_assert(source.size() == dimension_);
1387 : plumed_dbg_assert(sink.size() == dimension_);
1388 : // Start and end indices
1389 491 : index_t source_idx = getIndex(source);
1390 491 : index_t sink_idx = getIndex(sink);
1391 : // Path cost
1392 : double maximal_minimum = 0;
1393 : // In one dimension, path searching is very easy--either go one way if it's not periodic,
1394 : // or go both ways if it is periodic. There's no reason to pay the cost of Dijkstra.
1395 491 : if (dimension_ == 1) {
1396 : // Do a search from the grid source to grid sink that does not
1397 : // cross the grid boundary.
1398 147 : double curr_min_bias = getValue(source_idx);
1399 : // Either search from a high source to a low sink.
1400 147 : if (source_idx > sink_idx) {
1401 735 : for (index_t i = source_idx; i >= sink_idx; i--) {
1402 588 : if (curr_min_bias == 0.0) {
1403 : break;
1404 : }
1405 588 : curr_min_bias = fmin(curr_min_bias, getValue(i));
1406 : }
1407 : // Or search from a low source to a high sink.
1408 0 : } else if (source_idx < sink_idx) {
1409 0 : for (index_t i = source_idx; i <= sink_idx; i++) {
1410 0 : if (curr_min_bias == 0.0) {
1411 : break;
1412 : }
1413 0 : curr_min_bias = fmin(curr_min_bias, getValue(i));
1414 : }
1415 : }
1416 : maximal_minimum = curr_min_bias;
1417 : // If the grid is periodic, also do the search that crosses
1418 : // the grid boundary.
1419 147 : if (pbc_[0]) {
1420 42 : double c_min_bias = getValue(source_idx);
1421 : // Either go from a high source to the upper boundary and
1422 : // then from the bottom boundary to the sink
1423 42 : if (source_idx > sink_idx) {
1424 210 : for (index_t i = source_idx; i < maxsize_; i++) {
1425 168 : if (c_min_bias == 0.0) {
1426 : break;
1427 : }
1428 168 : c_min_bias = fmin(c_min_bias, getValue(i));
1429 : }
1430 210 : for (index_t i = 0; i <= sink_idx; i++) {
1431 168 : if (c_min_bias == 0.0) {
1432 : break;
1433 : }
1434 168 : c_min_bias = fmin(c_min_bias, getValue(i));
1435 : }
1436 : // Or go from a low source to the bottom boundary and
1437 : // then from the high boundary to the sink
1438 0 : } else if (source_idx < sink_idx) {
1439 0 : for (index_t i = source_idx; i > 0; i--) {
1440 0 : if (c_min_bias == 0.0) {
1441 : break;
1442 : }
1443 0 : c_min_bias = fmin(c_min_bias, getValue(i));
1444 : }
1445 0 : c_min_bias = fmin(c_min_bias, getValue(0));
1446 0 : for (index_t i = maxsize_ - 1; i <= sink_idx; i--) {
1447 0 : if (c_min_bias == 0.0) {
1448 : break;
1449 : }
1450 0 : c_min_bias = fmin(c_min_bias, getValue(i));
1451 : }
1452 : }
1453 : // If the boundary crossing paths was more biased, it's
1454 : // minimal bias replaces the non-boundary-crossing path's
1455 : // minimum.
1456 42 : maximal_minimum = fmax(maximal_minimum, c_min_bias);
1457 : }
1458 : // The one dimensional path search is complete.
1459 147 : return maximal_minimum;
1460 : // In two or more dimensions, path searching isn't trivial and we really
1461 : // do need to use a path search algorithm. Dijkstra is the simplest decent
1462 : // one. Using it we've never found the path search to be performance
1463 : // limiting in any solvated biomolecule test system, but faster options are
1464 : // easy to imagine if they become necessary. NB-In this case, we're actually
1465 : // using a greedy variant of Dijkstra's algorithm where the first possible
1466 : // path to a point always controls the path cost to that point. The structure
1467 : // of the cost function in this case guarantees that the calculated costs will
1468 : // be correct using this variant even though fine details of the paths may not
1469 : // match a normal Dijkstra search.
1470 344 : } else if (dimension_ > 1) {
1471 : // Prepare calculation temporaries for Dijkstra's algorithm.
1472 : // Minimal path costs from source to a given grid point
1473 344 : std::vector<double> mins_from_source = std::vector<double>(maxsize_, -1.0);
1474 : // Heap for tracking available steps, steps are recorded as std::pairs of
1475 : // an index and a value.
1476 344 : std::vector< std::pair<index_t, double> > next_steps;
1477 : std::pair<index_t, double> curr_indexed_val;
1478 344 : std::make_heap(next_steps.begin(), next_steps.end(), indexed_lt);
1479 : // The search begins at the source index.
1480 688 : next_steps.push_back(std::pair<index_t, double>(source_idx, getValue(source_idx)));
1481 344 : std::push_heap(next_steps.begin(), next_steps.end(), indexed_lt);
1482 : // At first no points have been examined and the optimal path has not been found.
1483 : bool path_not_found = true;
1484 : // Until a path is found,
1485 : while (path_not_found) {
1486 : // Examine the grid point currently most accessible from
1487 : // the set of all previously explored grid points by popping
1488 : // it from the top of the heap.
1489 29550 : std::pop_heap(next_steps.begin(), next_steps.end(), indexed_lt);
1490 : curr_indexed_val = next_steps.back();
1491 : next_steps.pop_back();
1492 : // Check if this point is the sink point, and if so
1493 : // finish the loop.
1494 29550 : if (curr_indexed_val.first == sink_idx) {
1495 : path_not_found = false;
1496 : maximal_minimum = curr_indexed_val.second;
1497 344 : break;
1498 : // Check if this point has reached the worst possible
1499 : // value, and if so stop looking for paths.
1500 29206 : } else if (curr_indexed_val.second == 0.0) {
1501 : maximal_minimum = 0.0;
1502 : break;
1503 : }
1504 : // If the search is not over, add this grid point's neighbors to the
1505 : // possible next points to search for the sink.
1506 29206 : std::vector<index_t> neighs = getNearestNeighbors(curr_indexed_val.first);
1507 144396 : for (unsigned k = 0; k < neighs.size(); k++) {
1508 115190 : index_t i = neighs[k];
1509 : // If the neighbor has not already been added to the list of possible next steps,
1510 115190 : if (mins_from_source[i] == -1.0) {
1511 : // Set the cost to reach it via a path through the current point being examined.
1512 43978 : mins_from_source[i] = fmin(curr_indexed_val.second, getValue(i));
1513 : // Add the neighboring point to the heap of potential next steps.
1514 43978 : next_steps.push_back(std::pair<index_t, double>(i, mins_from_source[i]));
1515 43978 : std::push_heap(next_steps.begin(), next_steps.end(), indexed_lt);
1516 : }
1517 : }
1518 : // Move on to the next best looking step along any of the paths
1519 : // growing from the source.
1520 : }
1521 : // The multidimensional path search is now complete.
1522 : return maximal_minimum;
1523 344 : }
1524 : return 0.0;
1525 : }
1526 :
1527 : }
|
