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 : #ifndef __PLUMED_tools_Grid_h
23 : #define __PLUMED_tools_Grid_h
24 :
25 : #include <vector>
26 : #include <string>
27 : #include <map>
28 : #include <cmath>
29 : #include <memory>
30 : #include <cstddef>
31 :
32 : namespace PLMD {
33 :
34 :
35 : // simple function to enable various weighting
36 :
37 : class WeightBase {
38 : public:
39 : virtual double projectInnerLoop(double &input, double &v)=0;
40 : virtual double projectOuterLoop(double &v)=0;
41 : virtual ~WeightBase() {}
42 : };
43 :
44 3 : class BiasWeight:public WeightBase {
45 : public:
46 : double beta,invbeta;
47 : double shift=0.0;
48 3 : explicit BiasWeight(double v) {
49 3 : beta=v;
50 3 : invbeta=1./beta;
51 : }
52 : //double projectInnerLoop(double &input, double &v) override {return input+exp(beta*v);}
53 : //double projectOuterLoop(double &v) override {return -invbeta*std::log(v);}
54 13603 : double projectInnerLoop(double &input, double &v) override {
55 13603 : auto betav=beta*v;
56 13603 : auto x=betav-shift;
57 13603 : if(x>0) {
58 1187 : shift=betav;
59 1187 : return input*std::exp(-x)+1;
60 : } else {
61 12416 : return input+std::exp(x);
62 : }
63 : }
64 187 : double projectOuterLoop(double &v) override {
65 187 : auto res=-invbeta*(std::log(v)+shift);
66 187 : shift=0.0;
67 187 : return res;
68 : }
69 : };
70 :
71 : class ProbWeight:public WeightBase {
72 : public:
73 : double beta,invbeta;
74 : explicit ProbWeight(double v) {
75 : beta=v;
76 : invbeta=1./beta;
77 : }
78 : double projectInnerLoop(double &input, double &v) override {
79 : return input+v;
80 : }
81 : double projectOuterLoop(double &v) override {
82 : return -invbeta*std::log(v);
83 : }
84 : };
85 :
86 :
87 :
88 :
89 :
90 :
91 : class Value;
92 : class IFile;
93 : class OFile;
94 : class KernelFunctions;
95 : class Communicator;
96 :
97 : /// \ingroup TOOLBOX
98 : class GridBase {
99 : public:
100 : // we use a size_t here
101 : // should be 8 bytes on all 64-bit machines
102 : // and more portable than "unsigned long long"
103 : typedef std::size_t index_t;
104 : // to restore old implementation (unsigned) use the following instead:
105 : // typedef unsigned index_t;
106 : /// Maximum dimension (exaggerated value).
107 : /// Can be used to replace local std::vectors with std::arrays (allocated on stack).
108 : static constexpr std::size_t maxdim=64;
109 : protected:
110 : std::string funcname;
111 : std::vector<std::string> argnames;
112 : std::vector<std::string> str_min_, str_max_;
113 : std::vector<double> min_,max_,dx_;
114 : std::vector<unsigned> nbin_;
115 : std::vector<bool> pbc_;
116 : index_t maxsize_;
117 : unsigned dimension_;
118 : bool dospline_, usederiv_;
119 : std::string fmt_; // format for output
120 : /// get "neighbors" for spline
121 : void getSplineNeighbors(const std::vector<unsigned> & indices, std::vector<index_t>& neigh, unsigned& nneigh )const;
122 : // std::vector<index_t> getSplineNeighbors(const std::vector<unsigned> & indices)const;
123 :
124 :
125 : public:
126 : /// this constructor here is Value-aware
127 : GridBase(const std::string& funcl, const std::vector<Value*> & args, const std::vector<std::string> & gmin,
128 : const std::vector<std::string> & gmax, const std::vector<unsigned> & nbin, bool dospline,
129 : bool usederiv);
130 : /// this constructor here is not Value-aware
131 : GridBase(const std::string& funcl, const std::vector<std::string> &names, const std::vector<std::string> & gmin,
132 : const std::vector<std::string> & gmax, const std::vector<unsigned> & nbin, bool dospline,
133 : bool usederiv, const std::vector<bool> &isperiodic, const std::vector<std::string> &pmin,
134 : const std::vector<std::string> &pmax );
135 : /// this is the real initializator
136 : void Init(const std::string & funcl, const std::vector<std::string> &names, const std::vector<std::string> & gmin,
137 : const std::vector<std::string> & gmax, const std::vector<unsigned> & nbin, bool dospline, bool usederiv,
138 : const std::vector<bool> &isperiodic, const std::vector<std::string> &pmin, const std::vector<std::string> &pmax);
139 : /// get lower boundary
140 : std::vector<std::string> getMin() const;
141 : /// get upper boundary
142 : std::vector<std::string> getMax() const;
143 : /// get bin size
144 : std::vector<double> getDx() const;
145 : double getDx(index_t j) const ;
146 : /// get bin volume
147 : double getBinVolume() const;
148 : /// get number of bins
149 : std::vector<unsigned> getNbin() const;
150 : /// get if periodic
151 : std::vector<bool> getIsPeriodic() const;
152 : /// get grid dimension
153 : unsigned getDimension() const;
154 : /// get argument names of this grid
155 : std::vector<std::string> getArgNames() const;
156 : /// get if the grid has derivatives
157 : bool hasDerivatives() const {
158 1986398 : return usederiv_;
159 : }
160 :
161 : /// methods to handle grid indices
162 : void getIndices(index_t index, std::vector<unsigned>& rindex) const;
163 : void getIndices(const std::vector<double> & x, std::vector<unsigned>& rindex) const;
164 : std::vector<unsigned> getIndices(index_t index) const;
165 : std::vector<unsigned> getIndices(const std::vector<double> & x) const;
166 : index_t getIndex(const std::vector<unsigned> & indices) const;
167 : index_t getIndex(const std::vector<double> & x) const;
168 : std::vector<double> getPoint(index_t index) const;
169 : std::vector<double> getPoint(const std::vector<unsigned> & indices) const;
170 : std::vector<double> getPoint(const std::vector<double> & x) const;
171 : /// faster versions relying on preallocated vectors
172 : void getPoint(index_t index,std::vector<double> & point) const;
173 : void getPoint(const std::vector<unsigned> & indices,std::vector<double> & point) const;
174 : void getPoint(const std::vector<double> & x,std::vector<double> & point) const;
175 :
176 : /// get neighbors
177 : std::vector<index_t> getNeighbors(index_t index,const std::vector<unsigned> & neigh) const;
178 : std::vector<index_t> getNeighbors(const std::vector<unsigned> & indices,const std::vector<unsigned> & neigh) const;
179 : std::vector<index_t> getNeighbors(const std::vector<double> & x,const std::vector<unsigned> & neigh) const;
180 : /// get nearest neighbors (those separated by exactly one lattice unit)
181 : std::vector<index_t> getNearestNeighbors(const index_t index) const;
182 : std::vector<index_t> getNearestNeighbors(const std::vector<unsigned> &indices) const;
183 :
184 : /// write header for grid file
185 : void writeHeader(OFile& file);
186 :
187 : /// read grid from file
188 : static std::unique_ptr<GridBase> create(const std::string&,const std::vector<Value*>&,IFile&,bool,bool,bool);
189 : /// read grid from file and check boundaries are what is expected from input
190 : static std::unique_ptr<GridBase> create(const std::string&,const std::vector<Value*>&, IFile&,
191 : const std::vector<std::string>&,const std::vector<std::string>&,
192 : const std::vector<unsigned>&,bool,bool,bool);
193 : /// get grid size
194 : virtual index_t getSize() const=0;
195 : /// get grid value
196 : virtual double getValue(index_t index) const=0;
197 : double getValue(const std::vector<unsigned> & indices) const;
198 : double getValue(const std::vector<double> & x) const;
199 : /// get grid value and derivatives
200 : virtual double getValueAndDerivatives(index_t index, std::vector<double>& der) const=0;
201 : double getValueAndDerivatives(const std::vector<unsigned> & indices, std::vector<double>& der) const;
202 : double getValueAndDerivatives(const std::vector<double> & x, std::vector<double>& der) const;
203 :
204 : /// set grid value
205 : virtual void setValue(index_t index, double value)=0;
206 : void setValue(const std::vector<unsigned> & indices, double value);
207 : /// set grid value and derivatives
208 : virtual void setValueAndDerivatives(index_t index, double value, std::vector<double>& der)=0;
209 : void setValueAndDerivatives(const std::vector<unsigned> & indices, double value, std::vector<double>& der);
210 : /// add to grid value
211 : virtual void addValue(index_t index, double value)=0;
212 : void addValue(const std::vector<unsigned> & indices, double value);
213 : /// add to grid value and derivatives
214 : virtual void addValueAndDerivatives(index_t index, double value, std::vector<double>& der)=0;
215 : void addValueAndDerivatives(const std::vector<unsigned> & indices, double value, std::vector<double>& der);
216 : /// add a kernel function to the grid
217 : void addKernel( const KernelFunctions& kernel );
218 :
219 : /// get minimum value
220 : virtual double getMinValue() const = 0;
221 : /// get maximum value
222 : virtual double getMaxValue() const = 0;
223 :
224 : /// dump grid on file
225 : virtual void writeToFile(OFile&)=0;
226 : /// dump grid to gaussian cube file
227 : void writeCubeFile(OFile&, const double& lunit);
228 :
229 3012 : virtual ~GridBase() {}
230 :
231 : /// set output format
232 : void setOutputFmt(const std::string & ss) {
233 179 : fmt_=ss;
234 197 : }
235 : /// reset output format to the default %14.9f format
236 : void resetToDefaultOutputFmt() {
237 : fmt_="%14.9f";
238 : }
239 : ///
240 : /// Find the maximum over paths of the minimum value of the gridded function along the paths
241 : /// for all paths of neighboring grid lattice points from a source point to a sink point.
242 : double findMaximalPathMinimum(const std::vector<double> &source, const std::vector<double> &sink);
243 : };
244 :
245 : class Grid : public GridBase {
246 : std::vector<double> grid_;
247 : std::vector<double> der_;
248 : double contour_location=0.0;
249 : public:
250 1323 : Grid(const std::string& funcl, const std::vector<Value*> & args, const std::vector<std::string> & gmin,
251 : const std::vector<std::string> & gmax,
252 1323 : const std::vector<unsigned> & nbin, bool dospline, bool usederiv):
253 1323 : GridBase(funcl,args,gmin,gmax,nbin,dospline,usederiv) {
254 1323 : grid_.assign(maxsize_,0.0);
255 1323 : if(usederiv_) {
256 191 : der_.assign(maxsize_*dimension_,0.0);
257 : }
258 1323 : }
259 : /// this constructor here is not Value-aware
260 128 : Grid(const std::string& funcl, const std::vector<std::string> &names, const std::vector<std::string> & gmin,
261 : const std::vector<std::string> & gmax, const std::vector<unsigned> & nbin, bool dospline,
262 : bool usederiv, const std::vector<bool> &isperiodic, const std::vector<std::string> &pmin,
263 128 : const std::vector<std::string> &pmax ):
264 128 : GridBase(funcl,names,gmin,gmax,nbin,dospline,usederiv,isperiodic,pmin,pmax) {
265 128 : grid_.assign(maxsize_,0.0);
266 128 : if(usederiv_) {
267 47 : der_.assign(maxsize_*dimension_,0.0);
268 : }
269 128 : }
270 : index_t getSize() const override;
271 : /// this is to access to Grid:: version of these methods (allowing overloading of virtual methods)
272 : using GridBase::getValue;
273 : using GridBase::getValueAndDerivatives;
274 : using GridBase::setValue;
275 : using GridBase::setValueAndDerivatives;
276 : using GridBase::addValue;
277 : using GridBase::addValueAndDerivatives;
278 : /// get grid value
279 : double getValue(index_t index) const override;
280 : /// get grid value and derivatives
281 : double getValueAndDerivatives(index_t index, std::vector<double>& der) const override;
282 :
283 : /// set grid value
284 : void setValue(index_t index, double value) override;
285 : /// set grid value and derivatives
286 : void setValueAndDerivatives(index_t index, double value, std::vector<double>& der) override;
287 : /// add to grid value
288 : void addValue(index_t index, double value) override;
289 : /// add to grid value and derivatives
290 : void addValueAndDerivatives(index_t index, double value, std::vector<double>& der) override;
291 :
292 : /// get minimum value
293 : double getMinValue() const override;
294 : /// get maximum value
295 : double getMaxValue() const override;
296 : /// Scale all grid values and derivatives by a constant factor
297 : void scaleAllValuesAndDerivatives( const double& scalef );
298 : /// Takes the scalef times the logarithm of all grid values and derivatives
299 : void logAllValuesAndDerivatives( const double& scalef );
300 : /// dump grid on file
301 : void writeToFile(OFile&) override;
302 :
303 : /// Set the minimum value of the grid to zero and translates accordingly
304 : void setMinToZero();
305 : /// apply function: takes pointer to function that accepts a double and apply
306 : void applyFunctionAllValuesAndDerivatives( double (*func)(double val), double (*funcder)(double valder) );
307 : /// Get the difference from the contour
308 : double getDifferenceFromContour(const std::vector<double> & x, std::vector<double>& der) const ;
309 : /// Find a set of points on a contour in the function
310 : void findSetOfPointsOnContour(const double& target, const std::vector<bool>& nosearch, unsigned& npoints, std::vector<std::vector<double> >& points );
311 : /// Since this method returns a concrete Grid, it should be here and not in GridBase - GB
312 : /// project a high dimensional grid onto a low dimensional one: this should be changed at some time
313 : /// to enable many types of weighting
314 : Grid project( const std::vector<std::string> & proj, WeightBase *ptr2obj );
315 : void projectOnLowDimension(double &val, std::vector<int> &varHigh, WeightBase* ptr2obj );
316 : void mpiSumValuesAndDerivatives( Communicator& comm );
317 : /// Integrate the function calculated on the grid
318 : double integrate( std::vector<unsigned>& npoints );
319 : void clear();
320 : };
321 :
322 :
323 : class SparseGrid : public GridBase {
324 :
325 : std::map<index_t,double> map_;
326 : std::map< index_t,std::vector<double> > der_;
327 :
328 : public:
329 6 : SparseGrid(const std::string& funcl, const std::vector<Value*> & args, const std::vector<std::string> & gmin,
330 : const std::vector<std::string> & gmax,
331 6 : const std::vector<unsigned> & nbin, bool dospline, bool usederiv):
332 6 : GridBase(funcl,args,gmin,gmax,nbin,dospline,usederiv) {}
333 :
334 : index_t getSize() const override;
335 : index_t getMaxSize() const;
336 :
337 : /// this is to access to Grid:: version of these methods (allowing overloading of virtual methods)
338 : using GridBase::getValue;
339 : using GridBase::getValueAndDerivatives;
340 : using GridBase::setValue;
341 : using GridBase::setValueAndDerivatives;
342 : using GridBase::addValue;
343 : using GridBase::addValueAndDerivatives;
344 :
345 : /// get grid value
346 : double getValue(index_t index) const override;
347 : /// get grid value and derivatives
348 : double getValueAndDerivatives(index_t index, std::vector<double>& der) const override;
349 :
350 : /// set grid value
351 : void setValue(index_t index, double value) override;
352 : /// set grid value and derivatives
353 : void setValueAndDerivatives(index_t index, double value, std::vector<double>& der) override;
354 : /// add to grid value
355 : void addValue(index_t index, double value) override;
356 : /// add to grid value and derivatives
357 : void addValueAndDerivatives(index_t index, double value, std::vector<double>& der) override;
358 :
359 : /// get minimum value
360 : double getMinValue() const override;
361 : /// get maximum value
362 : double getMaxValue() const override;
363 : /// dump grid on file
364 : void writeToFile(OFile&) override;
365 :
366 18 : virtual ~SparseGrid() {}
367 : };
368 : }
369 :
370 : #endif
|
