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