Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2018 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 :
23 : #include "Grid.h"
24 : #include "Tools.h"
25 : #include "core/Value.h"
26 : #include "File.h"
27 : #include "Exception.h"
28 : #include "KernelFunctions.h"
29 : #include "RootFindingBase.h"
30 : #include "Communicator.h"
31 :
32 : #include <vector>
33 : #include <cmath>
34 : #include <iostream>
35 : #include <sstream>
36 : #include <cstdio>
37 : #include <cfloat>
38 :
39 : using namespace std;
40 : namespace PLMD {
41 :
42 79 : Grid::Grid(const std::string& funcl, const std::vector<Value*> & args, const vector<std::string> & gmin,
43 79 : const vector<std::string> & gmax, const vector<unsigned> & nbin, bool dospline, bool usederiv, bool doclear) {
44 : // various checks
45 79 : plumed_massert(args.size()==gmin.size(),"grid min dimensions in input do not match number of arguments");
46 79 : plumed_massert(args.size()==nbin.size(),"number of bins on input do not match number of arguments");
47 79 : plumed_massert(args.size()==gmax.size(),"grid max dimensions in input do not match number of arguments");
48 79 : unsigned dim=gmax.size();
49 79 : std::vector<std::string> names;
50 158 : std::vector<bool> isperiodic;
51 158 : std::vector<string> pmin,pmax;
52 79 : names.resize( dim );
53 79 : isperiodic.resize( dim );
54 79 : pmin.resize( dim );
55 79 : pmax.resize( dim );
56 228 : for(unsigned int i=0; i<dim; ++i) {
57 149 : names[i]=args[i]->getName();
58 149 : if( args[i]->isPeriodic() ) {
59 37 : isperiodic[i]=true;
60 37 : args[i]->getDomain( pmin[i], pmax[i] );
61 : } else {
62 112 : isperiodic[i]=false;
63 112 : pmin[i]="0.";
64 112 : pmax[i]="0.";
65 : }
66 : }
67 : // this is a value-independent initializator
68 158 : Init(funcl,names,gmin,gmax,nbin,dospline,usederiv,doclear,isperiodic,pmin,pmax);
69 79 : }
70 :
71 2 : Grid::Grid(const std::string& funcl, const std::vector<string> &names, const std::vector<std::string> & gmin,
72 2 : const vector<std::string> & gmax, const std::vector<unsigned> & nbin, bool dospline, bool usederiv, bool doclear, const std::vector<bool> &isperiodic, const std::vector<std::string> &pmin, const std::vector<std::string> &pmax ) {
73 : // this calls the initializator
74 2 : Init(funcl,names,gmin,gmax,nbin,dospline,usederiv,doclear,isperiodic,pmin,pmax);
75 2 : }
76 :
77 81 : void Grid::Init(const std::string& funcl, const std::vector<std::string> &names, const vector<std::string> & gmin,
78 : const std::vector<std::string> & gmax, const std::vector<unsigned> & nbin, bool dospline, bool usederiv, bool doclear,
79 : const std::vector<bool> &isperiodic, const std::vector<std::string> &pmin, const std::vector<std::string> &pmax ) {
80 81 : contour_location=0.0; fmt_="%14.9f";
81 : // various checks
82 81 : plumed_massert(names.size()==gmin.size(),"grid dimensions in input do not match number of arguments");
83 81 : plumed_massert(names.size()==nbin.size(),"grid dimensions in input do not match number of arguments");
84 81 : plumed_massert(names.size()==gmax.size(),"grid dimensions in input do not match number of arguments");
85 81 : dimension_=gmax.size();
86 81 : str_min_=gmin; str_max_=gmax;
87 81 : argnames.resize( dimension_ );
88 81 : min_.resize( dimension_ );
89 81 : max_.resize( dimension_ );
90 81 : pbc_.resize( dimension_ );
91 232 : for(unsigned int i=0; i<dimension_; ++i) {
92 151 : argnames[i]=names[i];
93 151 : if( isperiodic[i] ) {
94 39 : pbc_[i]=true;
95 39 : str_min_[i]=pmin[i];
96 39 : str_max_[i]=pmax[i];
97 : } else {
98 112 : pbc_[i]=false;
99 : }
100 151 : Tools::convert(str_min_[i],min_[i]);
101 151 : Tools::convert(str_max_[i],max_[i]);
102 151 : funcname=funcl;
103 151 : plumed_massert(max_[i]>min_[i],"maximum in grid must be larger than minimum");
104 151 : plumed_massert(nbin[i]>0,"number of grid points must be greater than zero");
105 : }
106 81 : nbin_=nbin;
107 81 : dospline_=dospline;
108 81 : usederiv_=usederiv;
109 81 : if(dospline_) plumed_assert(dospline_==usederiv_);
110 81 : maxsize_=1;
111 232 : for(unsigned int i=0; i<dimension_; ++i) {
112 151 : dx_.push_back( (max_[i]-min_[i])/static_cast<double>( nbin_[i] ) );
113 151 : if( !pbc_[i] ) { max_[i] += dx_[i]; nbin_[i] += 1; }
114 151 : maxsize_*=nbin_[i];
115 : }
116 81 : if(doclear) clear();
117 81 : }
118 :
119 71 : void Grid::clear() {
120 71 : grid_.resize(maxsize_);
121 71 : if(usederiv_) der_.resize(maxsize_);
122 307017 : for(index_t i=0; i<maxsize_; ++i) {
123 306946 : grid_[i]=0.0;
124 306946 : if(usederiv_) {
125 306782 : (der_[i]).resize(dimension_);
126 306782 : for(unsigned int j=0; j<dimension_; ++j) der_[i][j]=0.0;
127 : }
128 : }
129 71 : }
130 :
131 276 : vector<std::string> Grid::getMin() const {
132 276 : return str_min_;
133 : }
134 :
135 276 : vector<std::string> Grid::getMax() const {
136 276 : return str_max_;
137 : }
138 :
139 7215276 : vector<double> Grid::getDx() const {
140 7215276 : return dx_;
141 : }
142 :
143 0 : double Grid::getBinVolume() const {
144 0 : double vol=1.;
145 0 : for(unsigned i=0; i<dx_.size(); ++i) vol*=dx_[i];
146 0 : return vol;
147 : }
148 :
149 13 : vector<bool> Grid::getIsPeriodic() const {
150 13 : return pbc_;
151 : }
152 :
153 13295 : vector<unsigned> Grid::getNbin() const {
154 13295 : return nbin_;
155 : }
156 :
157 176 : vector<string> Grid::getArgNames() const {
158 176 : return argnames;
159 : }
160 :
161 :
162 31066 : Grid::index_t Grid::getSize() const {
163 31066 : return maxsize_;
164 : }
165 :
166 5 : unsigned Grid::getDimension() const {
167 5 : return dimension_;
168 : }
169 :
170 : // we are flattening arrays using a column-major order
171 4841797 : Grid::index_t Grid::getIndex(const vector<unsigned> & indices) const {
172 : plumed_dbg_assert(indices.size()==dimension_);
173 14521525 : for(unsigned int i=0; i<dimension_; i++)
174 9679728 : if(indices[i]>=nbin_[i]) {
175 0 : std::string is;
176 0 : Tools::convert(i,is);
177 0 : std::string msg="ERROR: the system is looking for a value outside the grid along the " + is + " ("+getArgNames()[i]+")";
178 0 : plumed_merror(msg+" index!");
179 : }
180 4841797 : index_t index=indices[dimension_-1];
181 9679728 : for(unsigned int i=dimension_-1; i>0; --i) {
182 4837931 : index=index*nbin_[i-1]+indices[i-1];
183 : }
184 4841797 : return index;
185 : }
186 :
187 50284 : Grid::index_t Grid::getIndex(const vector<double> & x) const {
188 : plumed_dbg_assert(x.size()==dimension_);
189 50284 : return getIndex(getIndices(x));
190 : }
191 :
192 : // we are flattening arrays using a column-major order
193 4756053 : vector<unsigned> Grid::getIndices(index_t index) const {
194 4756053 : vector<unsigned> indices(dimension_);
195 4756053 : index_t kk=index;
196 4756053 : indices[0]=(index%nbin_[0]);
197 4756053 : for(unsigned int i=1; i<dimension_-1; ++i) {
198 0 : kk=(kk-indices[i-1])/nbin_[i-1];
199 0 : indices[i]=(kk%nbin_[i]);
200 : }
201 4756053 : if(dimension_>=2) {
202 4744558 : indices[dimension_-1]=((kk-indices[dimension_-2])/nbin_[dimension_-2]);
203 : }
204 4756053 : return indices;
205 : }
206 :
207 1858621 : vector<unsigned> Grid::getIndices(const vector<double> & x) const {
208 : plumed_dbg_assert(x.size()==dimension_);
209 1858621 : vector<unsigned> indices;
210 5573137 : for(unsigned int i=0; i<dimension_; ++i) {
211 3714516 : indices.push_back(unsigned(floor((x[i]-min_[i])/dx_[i])));
212 : }
213 1858621 : return indices;
214 : }
215 :
216 933007 : vector<double> Grid::getPoint(const vector<unsigned> & indices) const {
217 : plumed_dbg_assert(indices.size()==dimension_);
218 933007 : vector<double> x;
219 2790099 : for(unsigned int i=0; i<dimension_; ++i) {
220 1857092 : x.push_back(min_[i]+(double)(indices[i])*dx_[i]);
221 : }
222 933007 : return x;
223 : }
224 :
225 30636 : vector<double> Grid::getPoint(index_t index) const {
226 : plumed_dbg_assert(index<maxsize_);
227 30636 : return getPoint(getIndices(index));
228 : }
229 :
230 902371 : vector<double> Grid::getPoint(const vector<double> & x) const {
231 : plumed_dbg_assert(x.size()==dimension_);
232 902371 : return getPoint(getIndices(x));
233 : }
234 :
235 1095395 : void Grid::getPoint(index_t index,std::vector<double> & point) const {
236 : plumed_dbg_assert(index<maxsize_);
237 1095395 : getPoint(getIndices(index),point);
238 1095395 : }
239 :
240 1095395 : void Grid::getPoint(const std::vector<unsigned> & indices,std::vector<double> & point) const {
241 : plumed_dbg_assert(indices.size()==dimension_);
242 : plumed_dbg_assert(point.size()==dimension_);
243 3285018 : for(unsigned int i=0; i<dimension_; ++i) {
244 2189623 : point[i]=(min_[i]+(double)(indices[i])*dx_[i]);
245 : }
246 1095395 : }
247 :
248 0 : void Grid::getPoint(const std::vector<double> & x,std::vector<double> & point) const {
249 : plumed_dbg_assert(x.size()==dimension_);
250 0 : getPoint(getIndices(x),point);
251 0 : }
252 :
253 3595 : vector<Grid::index_t> Grid::getNeighbors
254 : (const vector<unsigned> &indices,const vector<unsigned> &nneigh)const {
255 : plumed_dbg_assert(indices.size()==dimension_ && nneigh.size()==dimension_);
256 :
257 3595 : vector<index_t> neighbors;
258 7190 : vector<unsigned> small_bin(dimension_);
259 :
260 3595 : unsigned small_nbin=1;
261 10764 : for(unsigned j=0; j<dimension_; ++j) {
262 7169 : small_bin[j]=(2*nneigh[j]+1);
263 7169 : small_nbin*=small_bin[j];
264 : }
265 :
266 7190 : vector<unsigned> small_indices(dimension_);
267 7190 : vector<unsigned> tmp_indices;
268 1175218 : for(unsigned index=0; index<small_nbin; ++index) {
269 1171623 : tmp_indices.resize(dimension_);
270 1171623 : unsigned kk=index;
271 1171623 : small_indices[0]=(index%small_bin[0]);
272 1171623 : for(unsigned i=1; i<dimension_-1; ++i) {
273 0 : kk=(kk-small_indices[i-1])/small_bin[i-1];
274 0 : small_indices[i]=(kk%small_bin[i]);
275 : }
276 1171623 : if(dimension_>=2) {
277 1170244 : small_indices[dimension_-1]=((kk-small_indices[dimension_-2])/small_bin[dimension_-2]);
278 : }
279 1171623 : unsigned ll=0;
280 3513490 : for(unsigned i=0; i<dimension_; ++i) {
281 2341867 : int i0=small_indices[i]-nneigh[i]+indices[i];
282 2341867 : if(!pbc_[i] && i0<0) continue;
283 2341707 : if(!pbc_[i] && i0>=nbin_[i]) continue;
284 2341655 : if( pbc_[i] && i0<0) i0=nbin_[i]-(-i0)%nbin_[i];
285 2341655 : if( pbc_[i] && i0>=nbin_[i]) i0%=nbin_[i];
286 2341655 : tmp_indices[ll]=((unsigned)i0);
287 2341655 : ll++;
288 : }
289 1171623 : tmp_indices.resize(ll);
290 1171623 : if(tmp_indices.size()==dimension_) {neighbors.push_back(getIndex(tmp_indices));}
291 : }
292 7190 : return neighbors;
293 : }
294 :
295 3595 : vector<Grid::index_t> Grid::getNeighbors
296 : (const vector<double> & x,const vector<unsigned> & nneigh)const {
297 : plumed_dbg_assert(x.size()==dimension_ && nneigh.size()==dimension_);
298 3595 : return getNeighbors(getIndices(x),nneigh);
299 : }
300 :
301 0 : vector<Grid::index_t> Grid::getNeighbors
302 : (index_t index,const vector<unsigned> & nneigh)const {
303 : plumed_dbg_assert(index<maxsize_ && nneigh.size()==dimension_);
304 0 : return getNeighbors(getIndices(index),nneigh);
305 : }
306 :
307 902371 : vector<Grid::index_t> Grid::getSplineNeighbors(const vector<unsigned> & indices)const {
308 : plumed_dbg_assert(indices.size()==dimension_);
309 902371 : vector<index_t> neighbors;
310 902371 : unsigned nneigh=unsigned(pow(2.0,int(dimension_)));
311 :
312 4509359 : for(unsigned int i=0; i<nneigh; ++i) {
313 3606988 : unsigned tmp=i;
314 3606988 : vector<unsigned> nindices;
315 10818468 : for(unsigned int j=0; j<dimension_; ++j) {
316 7211480 : unsigned i0=tmp%2+indices[j];
317 7211480 : tmp/=2;
318 7211480 : if(!pbc_[j] && i0==nbin_[j]) continue;
319 7211474 : if( pbc_[j] && i0==nbin_[j]) i0=0;
320 7211474 : nindices.push_back(i0);
321 : }
322 3606988 : if(nindices.size()==dimension_) neighbors.push_back(getIndex(nindices));
323 3606988 : }
324 902371 : return neighbors;
325 : }
326 :
327 0 : void Grid::addKernel( const KernelFunctions& kernel ) {
328 : plumed_dbg_assert( kernel.ndim()==dimension_ );
329 0 : std::vector<unsigned> nneighb=kernel.getSupport( dx_ );
330 0 : std::vector<index_t> neighbors=getNeighbors( kernel.getCenter(), nneighb );
331 0 : std::vector<double> xx( dimension_ ); std::vector<Value*> vv( dimension_ );
332 0 : std::string str_min, str_max;
333 0 : for(unsigned i=0; i<dimension_; ++i) {
334 0 : vv[i]=new Value();
335 0 : if( pbc_[i] ) {
336 0 : Tools::convert(min_[i],str_min);
337 0 : Tools::convert(max_[i],str_max);
338 0 : vv[i]->setDomain( str_min, str_max );
339 : } else {
340 0 : vv[i]->setNotPeriodic();
341 : }
342 : }
343 :
344 0 : std::vector<double> der( dimension_ );
345 0 : for(unsigned i=0; i<neighbors.size(); ++i) {
346 0 : index_t ineigh=neighbors[i];
347 0 : getPoint( ineigh, xx );
348 0 : for(unsigned j=0; j<dimension_; ++j) vv[j]->set(xx[j]);
349 0 : double newval = kernel.evaluate( vv, der, usederiv_ );
350 0 : if( usederiv_ ) addValueAndDerivatives( ineigh, newval, der );
351 0 : else addValue( ineigh, newval );
352 : }
353 :
354 0 : for(unsigned i=0; i<dimension_; ++i) delete vv[i];
355 0 : }
356 :
357 13457 : double Grid::getValue(index_t index) const {
358 : plumed_dbg_assert(index<maxsize_);
359 13457 : return grid_[index];
360 : }
361 :
362 0 : double Grid::getMinValue() const {
363 : double minval;
364 0 : minval=DBL_MAX;
365 0 : for(index_t i=0; i<grid_.size(); ++i) {
366 0 : if(grid_[i]<minval)minval=grid_[i];
367 : }
368 0 : return minval;
369 : }
370 :
371 18 : double Grid::getMaxValue() const {
372 : double maxval;
373 18 : maxval=DBL_MIN;
374 3618 : for(index_t i=0; i<grid_.size(); ++i) {
375 3600 : if(grid_[i]>maxval)maxval=grid_[i];
376 : }
377 18 : return maxval;
378 : }
379 :
380 :
381 13120 : double Grid::getValue(const vector<unsigned> & indices) const {
382 13120 : return getValue(getIndex(indices));
383 : }
384 :
385 1625 : double Grid::getValue(const vector<double> & x) const {
386 1625 : if(!dospline_) {
387 9 : return getValue(getIndex(x));
388 : } else {
389 1616 : vector<double> der(dimension_);
390 1616 : return getValueAndDerivatives(x,der);
391 : }
392 : }
393 :
394 3637074 : double Grid::getValueAndDerivatives
395 : (index_t index, vector<double>& der) const {
396 : plumed_dbg_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
397 3637074 : der=der_[index];
398 3637074 : return grid_[index];
399 : }
400 :
401 0 : double Grid::getValueAndDerivatives
402 : (const vector<unsigned> & indices, vector<double>& der) const {
403 0 : return getValueAndDerivatives(getIndex(indices),der);
404 : }
405 :
406 902371 : double Grid::getValueAndDerivatives
407 : (const vector<double> & x, vector<double>& der) const {
408 : plumed_dbg_assert(der.size()==dimension_ && usederiv_);
409 :
410 902371 : if(dospline_) {
411 : double X,X2,X3,value;
412 902371 : vector<double> fd(dimension_);
413 1804742 : vector<double> C(dimension_);
414 1804742 : vector<double> D(dimension_);
415 1804742 : vector<double> dder(dimension_);
416 : // reset
417 902371 : value=0.0;
418 902371 : for(unsigned int i=0; i<dimension_; ++i) der[i]=0.0;
419 :
420 1804742 : vector<unsigned> indices=getIndices(x);
421 1804742 : vector<index_t> neigh=getSplineNeighbors(indices);
422 1804742 : vector<double> xfloor=getPoint(x);
423 :
424 : // loop over neighbors
425 4509353 : for(unsigned int ipoint=0; ipoint<neigh.size(); ++ipoint) {
426 3606982 : double grid=getValueAndDerivatives(neigh[ipoint],dder);
427 3606982 : vector<unsigned> nindices=getIndices(neigh[ipoint]);
428 3606982 : double ff=1.0;
429 :
430 10818456 : for(unsigned j=0; j<dimension_; ++j) {
431 7211474 : int x0=1;
432 7211474 : if(nindices[j]==indices[j]) x0=0;
433 7211474 : double dx=getDx()[j];
434 7211474 : X=fabs((x[j]-xfloor[j])/dx-(double)x0);
435 7211474 : X2=X*X;
436 7211474 : X3=X2*X;
437 : double yy;
438 7211474 : if(fabs(grid)<0.0000001) yy=0.0;
439 1000192 : else yy=-dder[j]/grid;
440 7211474 : C[j]=(1.0-3.0*X2+2.0*X3) - (x0?-1.0:1.0)*yy*(X-2.0*X2+X3)*dx;
441 7211474 : D[j]=( -6.0*X +6.0*X2) - (x0?-1.0:1.0)*yy*(1.0-4.0*X +3.0*X2)*dx;
442 7211474 : D[j]*=(x0?-1.0:1.0)/dx;
443 7211474 : ff*=C[j];
444 : }
445 10818456 : for(unsigned j=0; j<dimension_; ++j) {
446 7211474 : fd[j]=D[j];
447 7211474 : for(unsigned i=0; i<dimension_; ++i) if(i!=j) fd[j]*=C[i];
448 : }
449 3606982 : value+=grid*ff;
450 3606982 : for(unsigned j=0; j<dimension_; ++j) der[j]+=grid*fd[j];
451 3606982 : }
452 1804742 : return value;
453 : } else {
454 0 : return getValueAndDerivatives(getIndex(x),der);
455 : }
456 : }
457 :
458 328 : void Grid::setValue(index_t index, double value) {
459 : plumed_dbg_assert(index<maxsize_ && !usederiv_);
460 328 : grid_[index]=value;
461 328 : }
462 :
463 0 : void Grid::setValue(const vector<unsigned> & indices, double value) {
464 0 : setValue(getIndex(indices),value);
465 0 : }
466 :
467 50275 : void Grid::setValueAndDerivatives
468 : (index_t index, double value, vector<double>& der) {
469 : plumed_dbg_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
470 50275 : grid_[index]=value;
471 50275 : der_[index]=der;
472 50275 : }
473 :
474 0 : void Grid::setValueAndDerivatives
475 : (const vector<unsigned> & indices, double value, vector<double>& der) {
476 0 : setValueAndDerivatives(getIndex(indices),value,der);
477 0 : }
478 :
479 0 : void Grid::addValue(index_t index, double value) {
480 : plumed_dbg_assert(index<maxsize_ && !usederiv_);
481 0 : grid_[index]+=value;
482 0 : }
483 :
484 0 : void Grid::addValue(const vector<unsigned> & indices, double value) {
485 0 : addValue(getIndex(indices),value);
486 0 : }
487 :
488 1151727 : void Grid::addValueAndDerivatives
489 : (index_t index, double value, vector<double>& der) {
490 : plumed_dbg_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
491 1151727 : grid_[index]+=value;
492 1151727 : for(unsigned int i=0; i<dimension_; ++i) der_[index][i]+=der[i];
493 1151727 : }
494 :
495 0 : void Grid::addValueAndDerivatives
496 : (const vector<unsigned> & indices, double value, vector<double>& der) {
497 0 : addValueAndDerivatives(getIndex(indices),value,der);
498 0 : }
499 :
500 6 : void Grid::scaleAllValuesAndDerivatives( const double& scalef ) {
501 6 : if(usederiv_) {
502 20968 : for(index_t i=0; i<grid_.size(); ++i) {
503 20962 : grid_[i]*=scalef;
504 20962 : for(unsigned j=0; j<dimension_; ++j) der_[i][j]*=scalef;
505 : }
506 : } else {
507 0 : for(index_t i=0; i<grid_.size(); ++i) grid_[i]*=scalef;
508 : }
509 6 : }
510 :
511 0 : void Grid::logAllValuesAndDerivatives( const double& scalef ) {
512 0 : if(usederiv_) {
513 0 : for(index_t i=0; i<grid_.size(); ++i) {
514 0 : grid_[i] = scalef*log(grid_[i]);
515 0 : for(unsigned j=0; j<dimension_; ++j) der_[i][j] = scalef/der_[i][j];
516 : }
517 : } else {
518 0 : for(index_t i=0; i<grid_.size(); ++i) grid_[i] = scalef*log(grid_[i]);
519 : }
520 0 : }
521 :
522 0 : void Grid::setMinToZero() {
523 0 : double min=grid_[0];
524 0 : for(index_t i=1; i<grid_.size(); ++i) if(grid_[i]<min) min=grid_[i];
525 0 : for(index_t i=0; i<grid_.size(); ++i) grid_[i] -= min;
526 0 : }
527 :
528 1 : void Grid::applyFunctionAllValuesAndDerivatives( double (*func)(double val), double (*funcder)(double valder) ) {
529 1 : if(usederiv_) {
530 901 : for(index_t i=0; i<grid_.size(); ++i) {
531 900 : grid_[i]=func(grid_[i]);
532 900 : for(unsigned j=0; j<dimension_; ++j) der_[i][j]=funcder(der_[i][j]);
533 : }
534 : } else {
535 0 : for(index_t i=0; i<grid_.size(); ++i) grid_[i]=func(grid_[i]);
536 : }
537 1 : }
538 :
539 98 : void Grid::writeHeader(OFile& ofile) {
540 282 : for(unsigned i=0; i<dimension_; ++i) {
541 184 : ofile.addConstantField("min_" + argnames[i]);
542 184 : ofile.addConstantField("max_" + argnames[i]);
543 184 : ofile.addConstantField("nbins_" + argnames[i]);
544 184 : ofile.addConstantField("periodic_" + argnames[i]);
545 : }
546 98 : }
547 :
548 98 : void Grid::writeToFile(OFile& ofile) {
549 98 : vector<double> xx(dimension_);
550 196 : vector<double> der(dimension_);
551 : double f;
552 98 : writeHeader(ofile);
553 30734 : for(index_t i=0; i<getSize(); ++i) {
554 30636 : xx=getPoint(i);
555 30636 : if(usederiv_) {f=getValueAndDerivatives(i,der);}
556 164 : else {f=getValue(i);}
557 30636 : if(i>0 && dimension_>1 && getIndices(i)[dimension_-2]==0) ofile.printf("\n");
558 84234 : for(unsigned j=0; j<dimension_; ++j) {
559 53598 : ofile.printField("min_" + argnames[j], str_min_[j] );
560 53598 : ofile.printField("max_" + argnames[j], str_max_[j] );
561 53598 : ofile.printField("nbins_" + argnames[j], static_cast<int>(nbin_[j]) );
562 53598 : if( pbc_[j] ) ofile.printField("periodic_" + argnames[j], "true" );
563 11510 : else ofile.printField("periodic_" + argnames[j], "false" );
564 : }
565 30636 : for(unsigned j=0; j<dimension_; ++j) { ofile.fmtField(" "+fmt_); ofile.printField(argnames[j],xx[j]); }
566 30636 : ofile.fmtField(" "+fmt_); ofile.printField(funcname,f);
567 30636 : if(usederiv_) for(unsigned j=0; j<dimension_; ++j) { ofile.fmtField(" "+fmt_); ofile.printField("der_" + argnames[j],der[j]); }
568 30636 : ofile.printField();
569 98 : }
570 98 : }
571 :
572 0 : void Grid::writeCubeFile(OFile& ofile, const double& lunit) {
573 0 : plumed_assert( dimension_==3 );
574 0 : ofile.printf("PLUMED CUBE FILE\n");
575 0 : ofile.printf("OUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z\n");
576 : // Number of atoms followed by position of origin (origin set so that center of grid is in center of cell)
577 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]));
578 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
579 0 : ofile.printf("%u %f %f %f\n",nbin_[1],0.0,lunit*dx_[1],0.0); // shape of voxel
580 0 : ofile.printf("%u %f %f %f\n",nbin_[2],0.0,0.0,lunit*dx_[2]);
581 0 : ofile.printf("%d %f %f %f\n",1,0.0,0.0,0.0); // Fake atom otherwise VMD doesn't work
582 0 : std::vector<unsigned> pp(3);
583 0 : for(pp[0]=0; pp[0]<nbin_[0]; ++pp[0]) {
584 0 : for(pp[1]=0; pp[1]<nbin_[1]; ++pp[1]) {
585 0 : for(pp[2]=0; pp[2]<nbin_[2]; ++pp[2]) {
586 0 : ofile.printf("%f ",getValue(pp) );
587 0 : if(pp[2]%6==5) ofile.printf("\n");
588 : }
589 0 : ofile.printf("\n");
590 : }
591 0 : }
592 0 : }
593 :
594 3 : Grid* Grid::create(const std::string& funcl, const std::vector<Value*> & args, IFile& ifile,
595 : const vector<std::string> & gmin,const vector<std::string> & gmax,
596 : const vector<unsigned> & nbin,bool dosparse, bool dospline, bool doder) {
597 3 : Grid* grid=Grid::create(funcl,args,ifile,dosparse,dospline,doder);
598 3 : std::vector<unsigned> cbin( grid->getNbin() );
599 6 : std::vector<std::string> cmin( grid->getMin() ), cmax( grid->getMax() );
600 9 : for(unsigned i=0; i<args.size(); ++i) {
601 6 : plumed_massert( cmin[i]==gmin[i], "mismatched grid min" );
602 6 : plumed_massert( cmax[i]==gmax[i], "mismatched grid max" );
603 6 : if( args[i]->isPeriodic() ) {
604 0 : plumed_massert( cbin[i]==nbin[i], "mismatched grid nbins" );
605 : } else {
606 6 : plumed_massert( (cbin[i]-1)==nbin[i], "mismatched grid nbins");
607 : }
608 : }
609 6 : return grid;
610 : }
611 :
612 5 : Grid* Grid::create(const std::string& funcl, const std::vector<Value*> & args, IFile& ifile, bool dosparse, bool dospline, bool doder)
613 : {
614 5 : Grid* grid=NULL;
615 5 : unsigned nvar=args.size(); bool hasder=false; std::string pstring;
616 10 : std::vector<int> gbin1(nvar); std::vector<unsigned> gbin(nvar);
617 10 : std::vector<std::string> labels(nvar),gmin(nvar),gmax(nvar);
618 10 : std::vector<std::string> fieldnames; ifile.scanFieldList( fieldnames );
619 : // Retrieve names for fields
620 5 : for(unsigned i=0; i<args.size(); ++i) labels[i]=args[i]->getName();
621 : // And read the stuff from the header
622 5 : plumed_massert( ifile.FieldExist( funcl ), "no column labelled " + funcl + " in in grid input");
623 14 : for(unsigned i=0; i<args.size(); ++i) {
624 9 : ifile.scanField( "min_" + labels[i], gmin[i]);
625 9 : ifile.scanField( "max_" + labels[i], gmax[i]);
626 9 : ifile.scanField( "periodic_" + labels[i], pstring );
627 9 : ifile.scanField( "nbins_" + labels[i], gbin1[i]);
628 9 : plumed_assert( gbin1[i]>0 );
629 9 : if( args[i]->isPeriodic() ) {
630 1 : plumed_massert( pstring=="true", "input value is periodic but grid is not");
631 2 : std::string pmin, pmax;
632 1 : args[i]->getDomain( pmin, pmax ); gbin[i]=gbin1[i];
633 2 : if( pmin!=gmin[i] || pmax!=gmax[i] ) plumed_merror("mismatch between grid boundaries and periods of values");
634 : } else {
635 8 : gbin[i]=gbin1[i]-1; // Note header in grid file indicates one more bin that there should be when data is not periodic
636 8 : plumed_massert( pstring=="false", "input value is not periodic but grid is");
637 : }
638 9 : hasder=ifile.FieldExist( "der_" + args[i]->getName() );
639 9 : if( doder && !hasder ) plumed_merror("missing derivatives from grid file");
640 13 : for(unsigned j=0; j<fieldnames.size(); ++j) {
641 17 : for(unsigned k=i+1; k<args.size(); ++k) {
642 4 : if( fieldnames[j]==labels[k] ) plumed_merror("arguments in input are not in same order as in grid file");
643 : }
644 13 : if( fieldnames[j]==labels[i] ) break;
645 : }
646 : }
647 :
648 5 : if(!dosparse) {grid=new Grid(funcl,args,gmin,gmax,gbin,dospline,doder);}
649 0 : else {grid=new SparseGrid(funcl,args,gmin,gmax,gbin,dospline,doder);}
650 :
651 10 : vector<double> xx(nvar),dder(nvar);
652 10 : vector<double> dx=grid->getDx();
653 : double f,x;
654 50285 : while( ifile.scanField(funcl,f) ) {
655 150625 : for(unsigned i=0; i<nvar; ++i) {
656 100350 : ifile.scanField(labels[i],x); xx[i]=x+dx[i]/2.0;
657 100350 : ifile.scanField( "min_" + labels[i], gmin[i]);
658 100350 : ifile.scanField( "max_" + labels[i], gmax[i]);
659 100350 : ifile.scanField( "nbins_" + labels[i], gbin1[i]);
660 100350 : ifile.scanField( "periodic_" + labels[i], pstring );
661 : }
662 50275 : if(hasder) { for(unsigned i=0; i<nvar; ++i) { ifile.scanField( "der_" + args[i]->getName(), dder[i] ); } }
663 50275 : index_t index=grid->getIndex(xx);
664 50275 : if(doder) {grid->setValueAndDerivatives(index,f,dder);}
665 0 : else {grid->setValue(index,f);}
666 50275 : ifile.scanField();
667 : }
668 10 : return grid;
669 : }
670 :
671 : // Sparse version of grid with map
672 0 : void SparseGrid::clear() {
673 0 : map_.clear();
674 0 : }
675 :
676 0 : Grid::index_t SparseGrid::getSize() const {
677 0 : return map_.size();
678 : }
679 :
680 0 : Grid::index_t SparseGrid::getMaxSize() const {
681 0 : return maxsize_;
682 : }
683 :
684 0 : double Grid::getDifferenceFromContour( const std::vector<double>& x, std::vector<double>& der ) const {
685 0 : return getValueAndDerivatives( x, der ) - contour_location;
686 : }
687 :
688 0 : void Grid::findSetOfPointsOnContour(const double& target, const std::vector<bool>& nosearch,
689 : unsigned& npoints, std::vector<std::vector<double> >& points ) {
690 : // Set contour location for function
691 0 : contour_location=target;
692 : // Resize points to maximum possible value
693 0 : points.resize( dimension_*maxsize_ );
694 :
695 : // Two points for search
696 0 : std::vector<unsigned> ind(dimension_);
697 0 : std::vector<double> direction( dimension_, 0 );
698 :
699 : // Run over whole grid
700 0 : npoints=0; RootFindingBase<Grid> mymin( this );
701 0 : for(unsigned i=0; i<maxsize_; ++i) {
702 0 : for(unsigned j=0; j<dimension_; ++j) ind[j]=getIndices(i)[j];
703 :
704 : // Get the value of a point on the grid
705 0 : double val1=getValue(i) - target;
706 :
707 : // Now search for contour in each direction
708 0 : bool edge=false;
709 0 : for(unsigned j=0; j<dimension_; ++j) {
710 0 : if( nosearch[j] ) continue ;
711 : // Make sure we don't search at the edge of the grid
712 0 : if( !pbc_[j] && (ind[j]+1)==nbin_[j] ) continue;
713 0 : else if( (ind[j]+1)==nbin_[j] ) { edge=true; ind[j]=0; }
714 0 : else ind[j]+=1;
715 0 : double val2=getValue(ind) - target;
716 0 : if( val1*val2<0 ) {
717 : // Use initial point location as first guess for search
718 0 : points[npoints].resize(dimension_); for(unsigned k=0; k<dimension_; ++k) points[npoints][k]=getPoint(i)[k];
719 : // Setup direction vector
720 0 : direction[j]=0.999999999*dx_[j];
721 : // And do proper search for contour point
722 0 : mymin.linesearch( direction, points[npoints], &Grid::getDifferenceFromContour );
723 0 : direction[j]=0.0; npoints++;
724 : }
725 0 : if( pbc_[j] && edge ) { edge=false; ind[j]=nbin_[j]-1; }
726 0 : else ind[j]-=1;
727 : }
728 0 : }
729 0 : }
730 :
731 0 : double SparseGrid::getValue(index_t index)const {
732 0 : plumed_assert(index<maxsize_);
733 0 : double value=0.0;
734 0 : iterator it=map_.find(index);
735 0 : if(it!=map_.end()) value=it->second;
736 0 : return value;
737 : }
738 :
739 380 : double SparseGrid::getValueAndDerivatives
740 : (index_t index, vector<double>& der)const {
741 380 : plumed_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
742 380 : double value=0.0;
743 380 : for(unsigned int i=0; i<dimension_; ++i) der[i]=0.0;
744 380 : iterator it=map_.find(index);
745 380 : if(it!=map_.end()) value=it->second;
746 380 : iterator_der itder=der_.find(index);
747 380 : if(itder!=der_.end()) der=itder->second;
748 380 : return value;
749 : }
750 :
751 0 : void SparseGrid::setValue(index_t index, double value) {
752 0 : plumed_assert(index<maxsize_ && !usederiv_);
753 0 : map_[index]=value;
754 0 : }
755 :
756 0 : void SparseGrid::setValueAndDerivatives
757 : (index_t index, double value, vector<double>& der) {
758 0 : plumed_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
759 0 : map_[index]=value;
760 0 : der_[index]=der;
761 0 : }
762 :
763 0 : void SparseGrid::addValue(index_t index, double value) {
764 0 : plumed_assert(index<maxsize_ && !usederiv_);
765 0 : map_[index]+=value;
766 0 : }
767 :
768 19684 : void SparseGrid::addValueAndDerivatives
769 : (index_t index, double value, vector<double>& der) {
770 19684 : plumed_assert(index<maxsize_ && usederiv_ && der.size()==dimension_);
771 19684 : map_[index]+=value;
772 19684 : der_[index].resize(dimension_);
773 19684 : for(unsigned int i=0; i<dimension_; ++i) der_[index][i]+=der[i];
774 19684 : }
775 :
776 0 : void SparseGrid::writeToFile(OFile& ofile) {
777 0 : vector<double> xx(dimension_);
778 0 : vector<double> der(dimension_);
779 : double f;
780 0 : writeHeader(ofile);
781 0 : ofile.fmtField(" "+fmt_);
782 0 : for(iterator it=map_.begin(); it!=map_.end(); ++it) {
783 0 : index_t i=(*it).first;
784 0 : xx=getPoint(i);
785 0 : if(usederiv_) {f=getValueAndDerivatives(i,der);}
786 0 : else {f=getValue(i);}
787 0 : if(i>0 && dimension_>1 && getIndices(i)[dimension_-2]==0) ofile.printf("\n");
788 0 : for(unsigned j=0; j<dimension_; ++j) {
789 0 : ofile.printField("min_" + argnames[j], str_min_[j] );
790 0 : ofile.printField("max_" + argnames[j], str_max_[j] );
791 0 : ofile.printField("nbins_" + argnames[j], static_cast<int>(nbin_[j]) );
792 0 : if( pbc_[j] ) ofile.printField("periodic_" + argnames[j], "true" );
793 0 : else ofile.printField("periodic_" + argnames[j], "false" );
794 : }
795 0 : for(unsigned j=0; j<dimension_; ++j) ofile.printField(argnames[j],xx[j]);
796 0 : ofile.printField(funcname, f);
797 0 : if(usederiv_) { for(unsigned j=0; j<dimension_; ++j) ofile.printField("der_" + argnames[j],der[j]); }
798 0 : ofile.printField();
799 0 : }
800 0 : }
801 :
802 :
803 13284 : void Grid::projectOnLowDimension(double &val, std::vector<int> &vHigh, WeightBase * ptr2obj ) {
804 13284 : unsigned i=0;
805 39688 : for(i=0; i<vHigh.size(); i++) {
806 26568 : if(vHigh[i]<0) { // this bin needs to be integrated out
807 : // parallelize here???
808 13284 : for(unsigned j=0; j<(getNbin())[i]; j++) {
809 13120 : vHigh[i]=int(j);
810 13120 : projectOnLowDimension(val,vHigh,ptr2obj); // recursive function: this is the core of the mechanism
811 13120 : vHigh[i]=-1;
812 : }
813 13448 : return; //
814 : }
815 : }
816 : // when there are no more bin to dig in then retrieve the value
817 13120 : if(i==vHigh.size()) {
818 : //std::cerr<<"POINT: ";
819 : //for(unsigned j=0;j<vHigh.size();j++){
820 : // std::cerr<<vHigh[j]<<" ";
821 : //}
822 13120 : std::vector<unsigned> vv(vHigh.size());
823 13120 : for(unsigned j=0; j<vHigh.size(); j++)vv[j]=unsigned(vHigh[j]);
824 : //
825 : // this is the real assignment !!!!! (hack this to have bias or other stuff)
826 : //
827 :
828 : // this case: produce fes
829 : //val+=exp(beta*getValue(vv)) ;
830 13120 : double myv=getValue(vv);
831 13120 : val=ptr2obj->projectInnerLoop(val,myv) ;
832 : // to be added: bias (same as before without negative sign)
833 : //std::cerr<<" VAL: "<<val <<endl;
834 : }
835 : }
836 :
837 2 : Grid Grid::project(const std::vector<std::string> & proj, WeightBase *ptr2obj ) {
838 : // find extrema only for the projection
839 4 : vector<string> smallMin,smallMax;
840 4 : vector<unsigned> smallBin;
841 4 : vector<unsigned> dimMapping;
842 4 : vector<bool> smallIsPeriodic;
843 4 : vector<string> smallName;
844 :
845 : // check if the two key methods are there
846 2 : WeightBase* pp = dynamic_cast<WeightBase*>(ptr2obj);
847 2 : if (!pp)plumed_merror("This WeightBase is not complete: you need a projectInnerLoop and projectOuterLoop ");
848 :
849 4 : for(unsigned j=0; j<proj.size(); j++) {
850 2 : for(unsigned i=0; i<getArgNames().size(); i++) {
851 2 : if(proj[j]==getArgNames()[i]) {
852 : unsigned offset;
853 : // note that at sizetime the non periodic dimension get a bin more
854 2 : if(getIsPeriodic()[i]) {offset=0;} else {offset=1;}
855 2 : smallMax.push_back(getMax()[i]);
856 2 : smallMin.push_back(getMin()[i]);
857 2 : smallBin.push_back(getNbin()[i]-offset);
858 2 : smallIsPeriodic.push_back(getIsPeriodic()[i]);
859 2 : dimMapping.push_back(i);
860 2 : smallName.push_back(getArgNames()[i]);
861 2 : break;
862 : }
863 : }
864 : }
865 2 : Grid smallgrid("projection",smallName,smallMin,smallMax,smallBin,false,false,true,smallIsPeriodic,smallMin,smallMax);
866 : // check that the two grids are commensurate
867 4 : for(unsigned i=0; i<dimMapping.size(); i++) {
868 2 : plumed_massert( (smallgrid.getMax())[i] == (getMax())[dimMapping[i]], "the two input grids are not compatible in max" );
869 2 : plumed_massert( (smallgrid.getMin())[i] == (getMin())[dimMapping[i]], "the two input grids are not compatible in min" );
870 2 : plumed_massert( (smallgrid.getNbin())[i]== (getNbin())[dimMapping[i]], "the two input grids are not compatible in bin" );
871 : }
872 4 : vector<unsigned> toBeIntegrated;
873 6 : for(unsigned i=0; i<getArgNames().size(); i++) {
874 4 : bool doappend=true;
875 6 : for(unsigned j=0; j<dimMapping.size(); j++) {
876 4 : if(dimMapping[j]==i) {doappend=false; break;}
877 : }
878 4 : if(doappend)toBeIntegrated.push_back(i);
879 : }
880 : //for(unsigned i=0;i<dimMapping.size();i++ ){
881 : // cerr<<"Dimension to preserve "<<dimMapping[i]<<endl;
882 : //}
883 : //for(unsigned i=0;i<toBeIntegrated.size();i++ ){
884 : // cerr<<"Dimension to integrate "<<toBeIntegrated[i]<<endl;
885 : //}
886 :
887 : // loop over all the points in the Grid, find the corresponding fixed index, rotate over all the other ones
888 166 : for(unsigned i=0; i<smallgrid.getSize(); i++) {
889 164 : std::vector<unsigned> v;
890 164 : v=smallgrid.getIndices(i);
891 328 : std::vector<int> vHigh((getArgNames()).size(),-1);
892 164 : for(unsigned j=0; j<dimMapping.size(); j++)vHigh[dimMapping[j]]=int(v[j]);
893 : // the vector vhigh now contains at the beginning the index of the low dimension and -1 for the values that need to be integrated
894 164 : double initval=0.;
895 164 : projectOnLowDimension(initval,vHigh, ptr2obj);
896 164 : smallgrid.setValue(i,initval);
897 164 : }
898 : // reset to zero just for biasing (this option can be evtl enabled in a future...)
899 : //double vmin;vmin=-smallgrid.getMinValue()+1;
900 166 : for(unsigned i=0; i<smallgrid.getSize(); i++) {
901 : // //if(dynamic_cast<BiasWeight*>(ptr2obj)){
902 : // // smallgrid.addValue(i,vmin);// go to 1
903 : // //}
904 164 : double vv=smallgrid.getValue(i);
905 164 : smallgrid.setValue(i,ptr2obj->projectOuterLoop(vv));
906 : // //if(dynamic_cast<BiasWeight*>(ptr2obj)){
907 : // // smallgrid.addValue(i,-vmin);// bring back to the value
908 : // //}
909 : }
910 :
911 4 : return smallgrid;
912 : }
913 :
914 0 : double Grid::integrate( std::vector<unsigned>& npoints ) {
915 0 : plumed_dbg_assert( npoints.size()==dimension_ ); plumed_assert( dospline_ );
916 :
917 0 : unsigned ntotgrid=1; double box_vol=1.0;
918 0 : std::vector<double> ispacing( npoints.size() );
919 0 : for(unsigned j=0; j<dimension_; ++j) {
920 0 : if( !pbc_[j] ) {
921 0 : ispacing[j] = ( max_[j] - dx_[j] - min_[j] ) / static_cast<double>( npoints[j] );
922 0 : npoints[j]+=1;
923 : } else {
924 0 : ispacing[j] = ( max_[j] - min_[j] ) / static_cast<double>( npoints[j] );
925 : }
926 0 : ntotgrid*=npoints[j]; box_vol*=ispacing[j];
927 : }
928 :
929 0 : std::vector<double> vals( dimension_ );
930 0 : std::vector<unsigned> t_index( dimension_ ); double integral=0.0;
931 0 : for(unsigned i=0; i<ntotgrid; ++i) {
932 0 : t_index[0]=(i%npoints[0]);
933 0 : unsigned kk=i;
934 0 : for(unsigned j=1; j<dimension_-1; ++j) { kk=(kk-t_index[j-1])/npoints[i-1]; t_index[j]=(kk%npoints[i]); }
935 0 : if( dimension_>=2 ) t_index[dimension_-1]=((kk-t_index[dimension_-1])/npoints[dimension_-2]);
936 :
937 0 : for(unsigned j=0; j<dimension_; ++j) vals[j]=min_[j] + t_index[j]*ispacing[j];
938 :
939 0 : integral += getValue( vals );
940 : }
941 :
942 0 : return box_vol*integral;
943 : }
944 :
945 0 : void Grid::mpiSumValuesAndDerivatives( Communicator& comm ) {
946 0 : comm.Sum( grid_ ); for(unsigned i=0; i<der_.size(); ++i) comm.Sum( der_[i] );
947 0 : }
948 :
949 2523 : }
|
