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 <functional>
27 : #include <array>
28 : #include <string>
29 : #include <map>
30 : #include <cmath>
31 : #include <memory>
32 : #include <cstddef>
33 :
34 : #include "Exception.h"
35 :
36 : namespace PLMD {
37 :
38 :
39 : // simple function to enable various weighting
40 :
41 : class WeightBase {
42 : public:
43 : virtual double projectInnerLoop(double &input, double &v)=0;
44 : virtual double projectOuterLoop(double &v)=0;
45 : virtual ~WeightBase() {}
46 : };
47 :
48 3 : class BiasWeight:public WeightBase {
49 : public:
50 : double beta,invbeta;
51 : double shift=0.0;
52 3 : explicit BiasWeight(double v) {
53 3 : beta=v;
54 3 : invbeta=1./beta;
55 : }
56 13603 : double projectInnerLoop(double &input, double &v) override {
57 13603 : auto betav=beta*v;
58 13603 : auto x=betav-shift;
59 13603 : if(x>0) {
60 1187 : shift=betav;
61 1187 : return input*std::exp(-x)+1;
62 : } else {
63 12416 : return input+std::exp(x);
64 : }
65 : }
66 187 : double projectOuterLoop(double &v) override {
67 187 : auto res=-invbeta*(std::log(v)+shift);
68 187 : shift=0.0;
69 187 : return res;
70 : }
71 : };
72 :
73 : class ProbWeight:public WeightBase {
74 : public:
75 : double beta,invbeta;
76 : explicit ProbWeight(double v) {
77 : beta=v;
78 : invbeta=1./beta;
79 : }
80 : double projectInnerLoop(double &input, double &v) override {
81 : return input+v;
82 : }
83 : double projectOuterLoop(double &v) override {
84 : return -invbeta*std::log(v);
85 : }
86 : };
87 :
88 :
89 :
90 :
91 :
92 :
93 : class Value;
94 : class IFile;
95 : class OFile;
96 : class KernelFunctions;
97 : class Communicator;
98 :
99 : /// \ingroup TOOLBOX
100 : class GridBase {
101 : public:
102 : // we use a size_t here
103 : // should be 8 bytes on all 64-bit machines
104 : // and more portable than "unsigned long long"
105 : typedef std::size_t index_t;
106 : // to restore old implementation (unsigned) use the following instead:
107 : // typedef unsigned index_t;
108 : /// Maximum dimension (exaggerated value).
109 : /// Can be used to replace local std::vectors with std::arrays (allocated on stack).
110 : static constexpr std::size_t maxdim=16;
111 :
112 :
113 : /**
114 : Auxiliary class for accelerating grid calculations.
115 :
116 : Many functions in the Grid class have loops running over the dimensions.
117 : These loops are always with a little number of iterations, so would likely
118 : benefit unrolling. Even better, if the compiler knew the dimension a priori
119 : it could remove the loops in most common cases (e.g. dimension=1 or 2).
120 :
121 : To obtain this, I added this AcceleratorBase class. Specifically:
122 : - Any method in the Grid class that is a bottleneck in performance and
123 : has a loop can be moved to this auxiliary class.
124 : - The method should be declared here as a pure virtual function
125 : of class AcceleratorBase. The implementation instead goes in
126 : class Accelerator
127 : - Class Accelerator is derived from AcceleratorBase but depends on
128 : dimension as a template variable.
129 : - The function AcceleratorBase::create takes case of constructing
130 : accelerators for each possible value of dimension.
131 :
132 : The trick works because we don't expect dimension to ever be very
133 : high. By restricting it to <=16, we can have a limited number of
134 : compiled instances. The compiler will choose at runtime which one should
135 : be used.
136 :
137 : This class, as well as the derived classes, are expected to be
138 : empty. In other words, they only contain the pointer to a virtual
139 : table that the compiler will use to call the proper implementation
140 : variant, with the correct dimension.
141 :
142 : \warning
143 : Interface might change at any time.
144 : Do not use this outside of class GridBase and children.
145 : */
146 : class AcceleratorBase {
147 : public:
148 : /// Creates an accelerator with proper dimension
149 : static std::unique_ptr<AcceleratorBase> create(unsigned dim);
150 : virtual ~AcceleratorBase() = default;
151 : /// Can be used to check which value of dimension was used
152 : virtual unsigned getDimension() const=0;
153 : virtual 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=0;
154 : virtual GridBase::index_t getIndex(const GridBase& grid, const std::vector<unsigned> & nbin_, const unsigned* indices,std::size_t indices_size) const=0;
155 : virtual 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=0;
156 : virtual void getIndices(const std::vector<unsigned> & nbin_, GridBase::index_t index, unsigned* indices, std::size_t indices_size) const=0;
157 : virtual 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=0;
158 : };
159 :
160 : /**
161 : Auxiliary class for managing AcceleratorBase.
162 :
163 : class GridBase should contain a std::unique_ptr<AcceleratorBase>, which
164 : can be used to access the specialized versions. However, this would
165 : make a GridBase object not copyable. Instead of redefining
166 : copy constructor and copy assignment for GridBase, which has a lot of
167 : members, we use this wrapper class that just manages the lifetime
168 : of the underlying Accelerator object.
169 :
170 : The underlying object is made accessible through the -> operator,
171 : so that one can simply call functions as accelerator->function.
172 :
173 : \warning
174 : Interface might change at any time.
175 : Do not use this outside of class GridBase and children.
176 : */
177 3242 : class AcceleratorHandler {
178 : /// This is the underlying pointer.
179 : std::unique_ptr<AcceleratorBase> ptr;
180 : public:
181 : /// Enable access to methods of AcceleratorBase
182 : const AcceleratorBase* operator->() const {
183 : return ptr.get();
184 : }
185 : /// Enable access to methods of AcceleratorBase
186 : AcceleratorBase* operator->() {
187 : return ptr.get();
188 : }
189 : /// Conversion to bool allows to check if the ptr has been set
190 : explicit operator bool() const noexcept {
191 : return bool(ptr);
192 : }
193 : /// Move constructor
194 66 : AcceleratorHandler(const AcceleratorHandler& other):
195 66 : ptr((other.ptr?AcceleratorBase::create(other->getDimension()):nullptr))
196 66 : {}
197 : /// Move assignment
198 1588 : AcceleratorHandler & operator=(const AcceleratorHandler & other) {
199 1588 : if(this!=&other) {
200 : ptr.reset();
201 1588 : if(other.ptr) {
202 3176 : ptr=AcceleratorBase::create(other->getDimension());
203 : }
204 : }
205 1588 : return *this;
206 : }
207 : /// Constructor without arguments result in a non-usable accelerator (dimension is unspecified)
208 : AcceleratorHandler() = default;
209 : /// Constructor with an argument creates an accelerator with a fixed dimensionality
210 1588 : AcceleratorHandler(unsigned dimension):
211 1588 : ptr(AcceleratorBase::create(dimension))
212 : {}
213 : };
214 :
215 : protected:
216 : AcceleratorHandler accelerator;
217 : std::string funcname;
218 : std::vector<std::string> argnames;
219 : std::vector<std::string> str_min_, str_max_;
220 : std::vector<double> min_,max_,dx_;
221 : std::vector<unsigned> nbin_;
222 : std::vector<bool> pbc_;
223 : index_t maxsize_;
224 : unsigned dimension_;
225 : bool dospline_, usederiv_;
226 : std::string fmt_; // format for output
227 : /// get "neighbors" for spline
228 : unsigned getSplineNeighbors(const unsigned* indices, std::size_t indices_size, index_t* neighbors, std::size_t neighbors_size)const;
229 : /// move the values after resize(): oldnbin are the numbers of points before resizing; the point with index k
230 : /// along direction i before resizing has index k+shift[i] after (it was removed if this is out of range)
231 : virtual void resizeData(const std::vector<unsigned> & oldnbin, const std::vector<long long> & shift);
232 :
233 :
234 : public:
235 : /// this constructor here is Value-aware
236 : GridBase(const std::string& funcl, const std::vector<Value*> & args, const std::vector<std::string> & gmin,
237 : const std::vector<std::string> & gmax, const std::vector<unsigned> & nbin, bool dospline,
238 : bool usederiv);
239 : /// this constructor here is not Value-aware
240 : GridBase(const std::string& funcl, const std::vector<std::string> &names, const std::vector<std::string> & gmin,
241 : const std::vector<std::string> & gmax, const std::vector<unsigned> & nbin, bool dospline,
242 : bool usederiv, const std::vector<bool> &isperiodic, const std::vector<std::string> &pmin,
243 : const std::vector<std::string> &pmax );
244 : /// this is the real initializator
245 : void Init(const std::string & funcl, const std::vector<std::string> &names, const std::vector<std::string> & gmin,
246 : const std::vector<std::string> & gmax, const std::vector<unsigned> & nbin, bool dospline, bool usederiv,
247 : const std::vector<bool> &isperiodic, const std::vector<std::string> &pmin, const std::vector<std::string> &pmax);
248 : /// get lower boundary
249 : std::vector<std::string> getMin() const;
250 : /// get upper boundary
251 : std::vector<std::string> getMax() const;
252 : /// get bin size
253 : std::vector<double> getDx() const;
254 : double getDx(index_t j) const ;
255 : /// get bin volume
256 : double getBinVolume() const;
257 : /// get number of bins
258 : std::vector<unsigned> getNbin() const;
259 : /// get if periodic
260 : std::vector<bool> getIsPeriodic() const;
261 : /// get grid dimension
262 : unsigned getDimension() const;
263 : /// true if x is within the range of the grid along all the non-periodic directions,
264 : /// i.e. if getIndices(x) gives a point of the grid
265 : bool isInside(const std::vector<double> & x) const;
266 : /// position of x along direction i, in units of the grid spacing, from the first point of the grid:
267 : /// its floor is the index of the point below x (also outside the grid, where it can be negative or too large)
268 : double getLatticeCoordinate(unsigned i, double x) const {
269 8490 : return (x-min_[i])/dx_[i];
270 : }
271 : /// Change the range of the grid along non-periodic directions: lower[i] (upper[i]) points are added at the lower (upper)
272 : /// end of direction i if positive, removed if negative. The spacing does not change, and the values (and derivatives)
273 : /// at the points that are kept do not change; new points are set to zero, the values at removed points are lost.
274 : /// The new bounds (getMin, getMax, written in the grid files) have 12 significant digits.
275 : void resize(const std::vector<long long> & lower, const std::vector<long long> & upper);
276 : /// get argument names of this grid
277 : std::vector<std::string> getArgNames() const;
278 : /// get if the grid has derivatives
279 : bool hasDerivatives() const {
280 1986449 : return usederiv_;
281 : }
282 :
283 : /// methods to handle grid indices
284 : void getIndices(index_t index, std::vector<unsigned>& rindex) const;
285 : void getIndices(index_t index, unsigned* rindex_data, std::size_t rindex_size) const;
286 : void getIndices(const std::vector<double> & x, std::vector<unsigned>& rindex) const;
287 : void getIndices(const std::vector<double> & x, unsigned* rindex_data,std::size_t rindex_size) const;
288 : std::vector<unsigned> getIndices(index_t index) const;
289 : std::vector<unsigned> getIndices(const std::vector<double> & x) const;
290 : index_t getIndex(const unsigned* indices,std::size_t indices_size) const;
291 : index_t getIndex(const std::vector<unsigned> & indices) const;
292 : index_t getIndex(const std::vector<double> & x) const;
293 : std::vector<double> getPoint(index_t index) const;
294 : std::vector<double> getPoint(const std::vector<unsigned> & indices) const;
295 : std::vector<double> getPoint(const std::vector<double> & x) const;
296 : /// faster versions relying on preallocated vectors
297 : void getPoint(index_t index,std::vector<double> & point) const;
298 : void getPoint(const std::vector<unsigned> & indices,std::vector<double> & point) const;
299 : void getPoint(const std::vector<double> & x,std::vector<double> & point) const;
300 : void getPoint(const unsigned* indices_data,std::size_t indices_size,std::vector<double> & point) const;
301 : void getPoint(const unsigned* indices_data,std::size_t indices_size,double* point,std::size_t point_size) const;
302 :
303 : /// get neighbors
304 : std::vector<index_t> getNeighbors(index_t index,const std::vector<unsigned> & neigh) const;
305 : std::vector<index_t> getNeighbors(const std::vector<unsigned> & indices,const std::vector<unsigned> & neigh) const;
306 : std::vector<index_t> getNeighbors(const std::vector<double> & x,const std::vector<unsigned> & neigh) const;
307 : /// get nearest neighbors (those separated by exactly one lattice unit)
308 : std::vector<index_t> getNearestNeighbors(const index_t index) const;
309 : std::vector<index_t> getNearestNeighbors(const std::vector<unsigned> &indices) const;
310 :
311 : /// write header for grid file
312 : void writeHeader(OFile& file);
313 :
314 : /// read grid from file
315 : static std::unique_ptr<GridBase> create(const std::string&,const std::vector<Value*>&,IFile&,bool,bool,bool);
316 : /// read grid from file and check boundaries are what is expected from input
317 : static std::unique_ptr<GridBase> create(const std::string&,const std::vector<Value*>&, IFile&,
318 : const std::vector<std::string>&,const std::vector<std::string>&,
319 : const std::vector<unsigned>&,bool,bool,bool);
320 : /// get grid size
321 : virtual index_t getSize() const=0;
322 : /// get grid value
323 : virtual double getValue(index_t index) const=0;
324 : double getValue(const std::vector<unsigned> & indices) const;
325 : double getValue(const std::vector<double> & x) const;
326 : /// get grid value and derivatives
327 : virtual double getValueAndDerivatives(index_t index, double* der, std::size_t der_size) const=0;
328 : double getValueAndDerivatives(index_t index, std::vector<double>& der) const;
329 : double getValueAndDerivatives(const std::vector<unsigned> & indices, std::vector<double>& der) const;
330 : double getValueAndDerivatives(const std::vector<double> & x, std::vector<double>& der) const;
331 :
332 : /// set grid value
333 : virtual void setValue(index_t index, double value)=0;
334 : void setValue(const std::vector<unsigned> & indices, double value);
335 : /// set grid value and derivatives
336 : virtual void setValueAndDerivatives(index_t index, double value, std::vector<double>& der)=0;
337 : void setValueAndDerivatives(const std::vector<unsigned> & indices, double value, std::vector<double>& der);
338 : /// add to grid value
339 : virtual void addValue(index_t index, double value)=0;
340 : void addValue(const std::vector<unsigned> & indices, double value);
341 : /// add to grid value and derivatives
342 : virtual void addValueAndDerivatives(index_t index, double value, std::vector<double>& der)=0;
343 : void addValueAndDerivatives(const std::vector<unsigned> & indices, double value, std::vector<double>& der);
344 : /// add to grid value and derivatives (der_size must be the dimension of the grid)
345 : virtual void addValueAndDerivatives(index_t index, double value, const double* der, std::size_t der_size);
346 : /// add a kernel function to the grid
347 : void addKernel( const KernelFunctions& kernel );
348 :
349 : /// call f(index,value) for each point that is stored, in increasing order of index: all the points of a dense grid,
350 : /// the points that have been set in a sparse grid (the others are zero)
351 : virtual void forEachStoredPoint(const std::function<void(index_t,double)> & f) const;
352 : /// get minimum value
353 : virtual double getMinValue() const = 0;
354 : /// get maximum value
355 : virtual double getMaxValue() const = 0;
356 :
357 : /// dump grid on file
358 : virtual void writeToFile(OFile&)=0;
359 : /// dump grid to gaussian cube file
360 : void writeCubeFile(OFile&, const double& lunit);
361 :
362 3308 : virtual ~GridBase() = default;
363 :
364 : /// set output format
365 : void setOutputFmt(const std::string & ss) {
366 730 : fmt_=ss;
367 266 : }
368 : /// reset output format to the default %14.9f format
369 : void resetToDefaultOutputFmt() {
370 : fmt_="%14.9f";
371 : }
372 : ///
373 : /// Find the maximum over paths of the minimum value of the gridded function along the paths
374 : /// for all paths of neighboring grid lattice points from a source point to a sink point.
375 : double findMaximalPathMinimum(const std::vector<double> &source, const std::vector<double> &sink);
376 : };
377 :
378 : class Grid : public GridBase {
379 : std::vector<double> grid_;
380 : std::vector<double> der_;
381 : double contour_location=0.0;
382 : public:
383 1445 : Grid(const std::string& funcl, const std::vector<Value*> & args, const std::vector<std::string> & gmin,
384 : const std::vector<std::string> & gmax,
385 1445 : const std::vector<unsigned> & nbin, bool dospline, bool usederiv):
386 1445 : GridBase(funcl,args,gmin,gmax,nbin,dospline,usederiv) {
387 1445 : grid_.assign(maxsize_,0.0);
388 1445 : if(usederiv_) {
389 313 : der_.assign(maxsize_*dimension_,0.0);
390 : }
391 1445 : }
392 : /// this constructor here is not Value-aware
393 128 : Grid(const std::string& funcl, const std::vector<std::string> &names, const std::vector<std::string> & gmin,
394 : const std::vector<std::string> & gmax, const std::vector<unsigned> & nbin, bool dospline,
395 : bool usederiv, const std::vector<bool> &isperiodic, const std::vector<std::string> &pmin,
396 128 : const std::vector<std::string> &pmax ):
397 128 : GridBase(funcl,names,gmin,gmax,nbin,dospline,usederiv,isperiodic,pmin,pmax) {
398 128 : grid_.assign(maxsize_,0.0);
399 128 : if(usederiv_) {
400 47 : der_.assign(maxsize_*dimension_,0.0);
401 : }
402 128 : }
403 : index_t getSize() const override;
404 : /// this is to access to Grid:: version of these methods (allowing overloading of virtual methods)
405 : using GridBase::getValue;
406 : using GridBase::getValueAndDerivatives;
407 : using GridBase::setValue;
408 : using GridBase::setValueAndDerivatives;
409 : using GridBase::addValue;
410 : using GridBase::addValueAndDerivatives;
411 : /// get grid value
412 : double getValue(index_t index) const override;
413 : /// get grid value and derivatives
414 : double getValueAndDerivatives(index_t index, double* der, std::size_t der_size) const override;
415 : /// set grid value
416 : void setValue(index_t index, double value) override;
417 : /// set grid value and derivatives
418 : void setValueAndDerivatives(index_t index, double value, std::vector<double>& der) override;
419 : /// add to grid value
420 : void addValue(index_t index, double value) override;
421 : /// add to grid value and derivatives
422 : void addValueAndDerivatives(index_t index, double value, std::vector<double>& der) override;
423 : void addValueAndDerivatives(index_t index, double value, const double* der, std::size_t der_size) override;
424 : protected:
425 : void resizeData(const std::vector<unsigned> & oldnbin, const std::vector<long long> & shift) override;
426 : public:
427 :
428 : /// get minimum value
429 : double getMinValue() const override;
430 : /// get maximum value
431 : double getMaxValue() const override;
432 : /// Scale all grid values and derivatives by a constant factor
433 : void scaleAllValuesAndDerivatives( const double& scalef );
434 : /// Takes the scalef times the logarithm of all grid values and derivatives
435 : void logAllValuesAndDerivatives( const double& scalef );
436 : /// dump grid on file
437 : void writeToFile(OFile&) override;
438 :
439 : /// Set the minimum value of the grid to zero and translates accordingly
440 : void setMinToZero();
441 : /// apply function: takes pointer to function that accepts a double and apply
442 : void applyFunctionAllValuesAndDerivatives( double (*func)(double val), double (*funcder)(double valder) );
443 : /// Get the difference from the contour
444 : double getDifferenceFromContour(const std::vector<double> & x, std::vector<double>& der) const ;
445 : /// Find a set of points on a contour in the function
446 : void findSetOfPointsOnContour(const double& target, const std::vector<bool>& nosearch, unsigned& npoints, std::vector<std::vector<double> >& points );
447 : /// Since this method returns a concrete Grid, it should be here and not in GridBase - GB
448 : /// project a high dimensional grid onto a low dimensional one: this should be changed at some time
449 : /// to enable many types of weighting
450 : Grid project( const std::vector<std::string> & proj, WeightBase *ptr2obj );
451 : void projectOnLowDimension(double &val, std::vector<int> &varHigh, WeightBase* ptr2obj );
452 : void mpiSumValuesAndDerivatives( Communicator& comm );
453 : /// Integrate the function calculated on the grid
454 : double integrate( std::vector<unsigned>& npoints );
455 : void clear();
456 : };
457 :
458 :
459 : class SparseGrid : public GridBase {
460 :
461 : std::map<index_t,double> map_;
462 : std::map< index_t,std::vector<double> > der_;
463 :
464 : public:
465 15 : SparseGrid(const std::string& funcl, const std::vector<Value*> & args, const std::vector<std::string> & gmin,
466 : const std::vector<std::string> & gmax,
467 15 : const std::vector<unsigned> & nbin, bool dospline, bool usederiv):
468 15 : GridBase(funcl,args,gmin,gmax,nbin,dospline,usederiv) {}
469 :
470 : index_t getSize() const override;
471 : index_t getMaxSize() const;
472 :
473 : /// this is to access to Grid:: version of these methods (allowing overloading of virtual methods)
474 : using GridBase::getValue;
475 : using GridBase::getValueAndDerivatives;
476 : using GridBase::setValue;
477 : using GridBase::setValueAndDerivatives;
478 : using GridBase::addValue;
479 : using GridBase::addValueAndDerivatives;
480 :
481 : /// get grid value
482 : double getValue(index_t index) const override;
483 : /// get grid value and derivatives
484 : double getValueAndDerivatives(index_t index, double* der, std::size_t der_size) const override;
485 :
486 : /// set grid value
487 : void setValue(index_t index, double value) override;
488 : /// set grid value and derivatives
489 : void setValueAndDerivatives(index_t index, double value, std::vector<double>& der) override;
490 : /// add to grid value
491 : void addValue(index_t index, double value) override;
492 : /// add to grid value and derivatives
493 : void addValueAndDerivatives(index_t index, double value, std::vector<double>& der) override;
494 : void addValueAndDerivatives(index_t index, double value, const double* der, std::size_t der_size) override;
495 : protected:
496 : void resizeData(const std::vector<unsigned> & oldnbin, const std::vector<long long> & shift) override;
497 : public:
498 : void forEachStoredPoint(const std::function<void(index_t,double)> & f) const override;
499 :
500 : /// get minimum value
501 : double getMinValue() const override;
502 : /// get maximum value
503 : double getMaxValue() const override;
504 : /// dump grid on file
505 : void writeToFile(OFile&) override;
506 :
507 45 : virtual ~SparseGrid() = default;
508 : };
509 :
510 :
511 : inline
512 6064176 : GridBase::index_t GridBase::getIndex(const unsigned* indices,std::size_t indices_size) const {
513 : plumed_dbg_assert(accelerator);
514 6064176 : return accelerator->getIndex(*this,nbin_,indices,indices_size);
515 : }
516 :
517 : inline
518 4791554 : GridBase::index_t GridBase::getIndex(const std::vector<unsigned> & indices) const {
519 : plumed_dbg_assert(indices.size()==dimension_);
520 4791554 : return getIndex(indices.data(),indices.size());
521 : }
522 :
523 : inline
524 1140094 : GridBase::index_t GridBase::getIndex(const std::vector<double> & x) const {
525 : plumed_dbg_assert(x.size()==dimension_);
526 : std::array<unsigned,maxdim> indices;
527 1140094 : getIndices(x,indices.data(),dimension_);
528 1140094 : return getIndex(indices.data(),dimension_);
529 : }
530 :
531 : inline
532 7904699 : std::vector<unsigned> GridBase::getIndices(index_t index) const {
533 7904699 : std::vector<unsigned> indices(dimension_);
534 7904699 : getIndices(index,indices.data(),indices.size());
535 7904699 : return indices;
536 : }
537 :
538 : inline
539 37954913 : void GridBase::getIndices(index_t index, unsigned* indices, std::size_t indices_size) const {
540 37954913 : plumed_assert(indices_size==dimension_);
541 37954913 : plumed_assert(accelerator);
542 37954913 : accelerator->getIndices(nbin_,index,indices,dimension_);
543 37954913 : }
544 :
545 : inline
546 : void GridBase::getIndices(index_t index, std::vector<unsigned>& indices) const {
547 : if (indices.size()!=dimension_) {
548 : indices.resize(dimension_);
549 : }
550 : getIndices(index,indices.data(),indices.size());
551 : }
552 :
553 : inline
554 : std::vector<unsigned> GridBase::getIndices(const std::vector<double> & x) const {
555 : plumed_dbg_assert(x.size()==dimension_);
556 : std::vector<unsigned> indices(dimension_);
557 : getIndices(x,indices);
558 : return indices;
559 : }
560 :
561 : inline
562 11692 : void GridBase::getIndices(const std::vector<double> & x, std::vector<unsigned>& indices) const {
563 11692 : indices.resize(dimension_);
564 11692 : getIndices(x,indices.data(),indices.size());
565 11692 : }
566 :
567 : inline
568 1200058 : void GridBase::getIndices(const std::vector<double> & x, unsigned* rindex_data,std::size_t rindex_size) const {
569 : plumed_dbg_assert(x.size()==dimension_);
570 : plumed_dbg_assert(rindex_size==dimension_);
571 1200058 : accelerator->getIndices(min_,dx_,x,rindex_data,rindex_size);
572 1200058 : }
573 :
574 : inline
575 33313933 : void GridBase::getPoint(const unsigned* indices,std::size_t indices_size,double* point,std::size_t point_size) const {
576 : plumed_dbg_assert(indices_size==dimension_);
577 : plumed_dbg_assert(point_size==dimension_);
578 : plumed_dbg_assert(accelerator);
579 33313933 : accelerator->getPoint(min_,dx_,indices,indices_size,point,point_size);
580 33313933 : }
581 :
582 : inline
583 2720117 : std::vector<double> GridBase::getPoint(const std::vector<unsigned> & indices) const {
584 : plumed_dbg_assert(indices.size()==dimension_);
585 2720117 : std::vector<double> x(dimension_);
586 2720117 : getPoint(indices,x);
587 2720117 : return x;
588 : }
589 :
590 : inline
591 28851721 : std::vector<double> GridBase::getPoint(index_t index) const {
592 : plumed_dbg_assert(index<maxsize_);
593 28851721 : std::vector<double> x(dimension_);
594 28851721 : getPoint(index,x);
595 28851721 : return x;
596 : }
597 :
598 : inline
599 : std::vector<double> GridBase::getPoint(const std::vector<double> & x) const {
600 : plumed_dbg_assert(x.size()==dimension_);
601 : std::vector<double> r(dimension_);
602 : getPoint(x,r);
603 : return r;
604 : }
605 :
606 : inline
607 29917686 : void GridBase::getPoint(index_t index,std::vector<double> & point) const {
608 : plumed_dbg_assert(index<maxsize_);
609 : std::array<unsigned,maxdim> indices;
610 29917686 : getIndices(index,indices.data(),dimension_);
611 29917686 : getPoint(indices.data(),dimension_,point);
612 29917686 : }
613 :
614 : inline
615 2720117 : void GridBase::getPoint(const std::vector<unsigned> & indices,std::vector<double> & point) const {
616 : plumed_dbg_assert(indices.size()==dimension_);
617 : plumed_dbg_assert(point.size()==dimension_);
618 2720117 : getPoint(indices.data(),indices.size(),point.data(),point.size());
619 2720117 : }
620 :
621 : inline
622 29917686 : void GridBase::getPoint(const unsigned* indices_data,std::size_t indices_size,std::vector<double> & point) const {
623 : plumed_dbg_assert(indices_size==dimension_);
624 : plumed_dbg_assert(point.size()==dimension_);
625 29917686 : getPoint(indices_data,indices_size,point.data(),point.size());
626 29917686 : }
627 :
628 : inline
629 : void GridBase::getPoint(const std::vector<double> & x,std::vector<double> & point) const {
630 : plumed_dbg_assert(x.size()==dimension_);
631 : std::array<unsigned,maxdim> indices;
632 : getIndices(x,indices.data(),dimension_);
633 : getPoint(indices.data(),dimension_,point.data(),point.size());
634 : }
635 :
636 : inline
637 : std::vector<GridBase::index_t> GridBase::getNeighbors(const std::vector<unsigned> &indices,const std::vector<unsigned> &nneigh)const {
638 : plumed_dbg_assert(accelerator);
639 : return accelerator->getNeighbors(*this,nbin_,pbc_,indices.data(),indices.size(),nneigh);
640 : }
641 :
642 : inline
643 3451 : std::vector<GridBase::index_t> GridBase::getNeighbors(const std::vector<double> & x,const std::vector<unsigned> & nneigh)const {
644 : plumed_dbg_assert(x.size()==dimension_ && nneigh.size()==dimension_);
645 : std::array<unsigned,maxdim> indices;
646 : plumed_dbg_assert(accelerator);
647 3451 : accelerator->getIndices(min_,dx_,x,indices.data(),dimension_);
648 3451 : return accelerator->getNeighbors(*this,nbin_,pbc_,indices.data(),dimension_,nneigh);
649 : }
650 :
651 : inline
652 58413 : std::vector<GridBase::index_t> GridBase::getNeighbors(index_t index,const std::vector<unsigned> & nneigh)const {
653 : plumed_dbg_assert(index<maxsize_ && nneigh.size()==dimension_);
654 : std::array<unsigned,maxdim> indices;
655 : plumed_dbg_assert(accelerator);
656 58413 : accelerator->getIndices(nbin_,index,indices.data(),dimension_);
657 58413 : return accelerator->getNeighbors(*this,nbin_,pbc_,indices.data(),dimension_,nneigh);
658 : }
659 :
660 :
661 :
662 :
663 :
664 : }
665 :
666 : #endif
|
