SSAGES  0.8.3
Software Suite for Advanced General Ensemble Simulations
GridBase.h
1 
20 #pragma once
21 
22 #include <cmath>
23 #include <exception>
24 #include <vector>
25 
26 namespace SSAGES
27 {
28 
30 
37 template<typename T>
38 class GridBase
39 {
40 protected:
42  std::vector<T> data_;
43 
45  size_t dimension_;
46 
48  std::vector<int> numPoints_;
49 
51  std::pair< std::vector<double>,std::vector<double> > edges_;
52 
54  std::vector<bool> isPeriodic_;
55 
57  std::vector<int> wrapIndices(const std::vector<int> &indices) const
58  {
59  std::vector<int> newIndices(indices);
60  for (size_t i = 0; i < dimension_; ++i) {
61  if (!GetPeriodic(i)) {
62  continue;
63  }
64 
65  int index = indices.at(i);
66  while (index < 0) { index += numPoints_[i]; }
67  while (index >= numPoints_[i]) { index -= numPoints_[i]; }
68  if(index == numPoints_[i]-1)
69  index = 0;
70  newIndices.at(i) = index;
71  }
72 
73  return newIndices;
74  }
75 
77 
85  virtual size_t mapTo1d(const std::vector<int> &indices) const = 0;
86 
87 protected:
89 
99  GridBase(std::vector<int> numPoints,
100  std::vector<double> lower,
101  std::vector<double> upper,
102  std::vector<bool> isPeriodic)
103  : dimension_(numPoints.size()),
104  numPoints_(numPoints),
105  edges_(std::pair< std::vector<double>, std::vector<double> >(lower, upper)),
106  isPeriodic_(isPeriodic)
107  {
108  // Check that vector sizes are correct
109  if (edges_.first.size() != dimension_ ||
110  edges_.second.size() != dimension_) {
111  throw std::invalid_argument("Size of vector containing upper or "
112  "lower edges, does not match size of vector containing "
113  "number of grid points.");
114  }
115  if (isPeriodic_.size() == 0) {
116  // Default: Non-periodic in all dimensions
117  isPeriodic.resize(dimension_, false);
118  } else if (isPeriodic_.size() != dimension_) {
119  throw std::invalid_argument("Size of vector isPeriodic does not "
120  "match size of vector containing number of grid points.");
121  }
122  for(size_t i=0 ; i < isPeriodic_.size() ; i++)
123  {
124  if(isPeriodic_[i])
125  {
126  if(numPoints_[i] <= 1)
127  {
128  throw std::invalid_argument("A periodic grid is incompatible with a grid size of 1.");
129  }
130  double spacing = (edges_.second[i] - edges_.first[i]) / (numPoints_[i]-1);
131  edges_.first[i] -= spacing/2;
132  edges_.second[i] += spacing/2;
133  }
134  }
135 
136  }
137 
138 public:
139 
141  void syncGrid()
142  {
143  //Convenience
144  size_t dim = this->GetDimension();
145 
146  //Preallocate
147  std::vector<int> navigate(dim);
148 
149  //Loop over all surfaces. Number of surfaces = dim.
150  for(size_t i = 0 ; i < dim ; i++)
151  {
152  //Check if periodic in this dimension.
153  if(isPeriodic_[i])
154  {
155 
156  //Calculate surface size. This is equal to number of points in each other dimension multiplied.
157  size_t surfsize = 1;
158  for(size_t j = 0 ; j < dim ; j++)
159  {
160  if(i != j)
161  surfsize*=numPoints_[j];
162  }
163 
164  //Loop over all points of this surface on the 0 side and copy to end.
165  for(size_t j = 0 ; j < surfsize ; j++)
166  {
167  int runningcount = j;
168  for(size_t k = 0 ; k < dim ; k++)
169  {
170  if(k == i)
171  navigate[k] = 0;
172  else
173  {
174  navigate[k] = runningcount%numPoints_[k];
175  runningcount = runningcount / numPoints_[k];
176  }
177  }
178  auto temp = this->at(navigate);
179  navigate[i] = numPoints_[i]-1;
180  this->at(navigate) = temp;
181  }
182  }
183  }
184  }
185 
187 
190  size_t GetDimension() const
191  {
192  return dimension_;
193  }
194 
196 
200  const std::vector<int> GetNumPoints() const
201  {
202  return numPoints_;
203  }
204 
206 
212  int GetNumPoints(size_t dim) const
213  {
214  if (dim >= GetDimension()) {
215  std::cerr << "Warning! Grid size requested for a dimension larger "
216  "than the grid dimensionality!\n";
217  return 0;
218  }
219  return numPoints_.at(dim);
220  }
221 
223 
226  const std::vector<double> GetLower() const
227  {
228  std::vector<double> lower(dimension_);
229  for(size_t i = 0; i<dimension_;++i)
230  if(GetPeriodic(i)) {lower[i] = edges_.first[i] - ((edges_.first[i] - edges_.second[i]) / (numPoints_[i]))/2;}
231  else {lower[i] = edges_.first[i];}
232  return lower;
233  }
234 
236 
242  double GetLower(size_t dim) const
243  {
244  if (dim >= GetDimension()) {
245  std::cerr << "Warning! Lower edge requested for a dimension larger "
246  "than the grid dimensionality!\n";
247  return 0.0;
248  }
249  if(GetPeriodic(dim)){return edges_.first[dim] - ((edges_.first[dim] - edges_.second[dim]) / (numPoints_[dim]))/2;}
250  else{return edges_.first[dim];}
251  }
252 
254 
257  const std::vector<double> GetUpper() const
258  {
259  std::vector<double> upper(dimension_);
260  for(size_t i = 0; i<dimension_;++i)
261  if(GetPeriodic(i)) {upper[i] = edges_.second[i] + ((edges_.first[i] - edges_.second[i]) / (numPoints_[i]))/2;}
262  else {upper[i] = edges_.second[i];}
263  return upper;
264  }
265 
266 
268 
274  double GetUpper(size_t dim) const
275  {
276  if (dim >= GetDimension()) {
277  std::cerr << "Warning! Upper edge requested for a dimension larger "
278  "than the grid dimensionality!\n";
279  return 0.0;
280  }
281  if(GetPeriodic(dim)){return edges_.second[dim] + ((edges_.first[dim] - edges_.second[dim]) / (numPoints_[dim]))/2;}
282  else{return edges_.second[dim];}
283  }
284 
286 
290  const std::vector<bool>& GetPeriodic() const
291  {
292  return isPeriodic_;
293  }
294 
296 
303  bool GetPeriodic(size_t dim) const
304  {
305  if (dim >= GetDimension()) {
306  std::cerr << "Warning! Periodicity requested for a dimension larger "
307  "than the grid dimensionality!\n";
308  return false;
309  }
310  return GetPeriodic().at(dim);
311  }
312 
314 
321  size_t size() const
322  {
323  return data_.size();
324  }
325 
327 
333  T *data()
334  {
335  return data_.data();
336  }
337 
339 
345  T const *data() const
346  {
347  return data_.data();
348  }
349 
351 
367  std::vector<int> GetIndices(const std::vector<double> &x) const
368  {
369  // Check that input vector has the correct dimensionality
370  if (x.size() != dimension_) {
371  throw std::invalid_argument("Specified point has a larger "
372  "dimensionality than the grid.");
373  }
374 
375  std::vector<int> indices(dimension_);
376  for (size_t i = 0; i < dimension_; ++i)
377  {
378  double xpos = x.at(i);
379  if (!GetPeriodic(i))
380  {
381  if (xpos < edges_.first[i])
382  {
383  indices.at(i) = -1;
384  continue;
385  }
386  else if (xpos > edges_.second[i])
387  {
388  indices.at(i) = numPoints_[i];
389  continue;
390  }
391 
392  }
393 
394  // To make sure, the value is rounded in the correct direction.
395  double spacing = (edges_.second[i] - edges_.first[i]) / (numPoints_[i]);
396  indices.at(i) = std::floor( (xpos - edges_.first[i]) / spacing);
397  }
398 
399  return wrapIndices(indices);
400  }
401 
403 
458 
467  int GetIndex(double x) const
468  {
469  if (dimension_ != 1) {
470  throw std::invalid_argument("1d Grid index can only be accessed for "
471  "1d-Grids can be accessed with a.");
472  }
473  return GetIndices({x}).at(0);
474  }
475 
476 
478  /*/!
479  *\param coords
480  *\param dim
481  * To be finalized.
482 
483  std::vector<double> GetDist(const std::vector<double> &coords)
484  {
485 
486  }*/
487 
488 
490 
498  std::vector<double> GetCoordinates(const std::vector<int> &indices)
499  {
500  if (indices.size() != dimension_) {
501  throw std::invalid_argument(
502  "Grid indices specified for GetCoordinates() do not have "
503  "the same dimensionality as the grid.");
504  }
505 
506  std::vector<double> v(dimension_);
507 
508  for (size_t i = 0; i < dimension_; ++i) {
509  double spacing = (edges_.second[i] - edges_.first[i]) / (numPoints_[i]);
510  v.at(i) = edges_.first[i] + (indices[i] + 0.5)*spacing;
511  }
512 
513  return v;
514  }
515 
517 
523  double GetCoordinate(int index)
524  {
525  return GetCoordinates({index}).at(0);
526  }
527 
529 
541  const T& at(const std::vector<int> &indices) const
542  {
543  // Check that indices are in bound.
544  if (indices.size() != GetDimension()) {
545  throw std::invalid_argument("Dimension of indices does not match "
546  "dimension of the grid.");
547  }
548 
549  return data_.at(mapTo1d(indices));
550  }
551 
553 
557  T& at(const std::vector<int> &indices)
558  {
559  return const_cast<T&>(static_cast<const GridBase<T>* >(this)->at(indices));
560  }
561 
563 
573  template<typename R>
574  const T& at(std::initializer_list<R>&& x) const
575  {
576  return at(static_cast<std::vector<R> >(x));
577  }
578 
580 
590  template<typename R>
591  T& at(std::initializer_list<R>&& x)
592  {
593  return at(static_cast<std::vector<R> >(x));
594  }
595 
597 
603  const T& at(int index) const
604  {
605  if (dimension_ != 1) {
606  throw std::invalid_argument("Only 1d-Grids can be accessed with a "
607  "single integer as the index.");
608  }
609  return at({index});
610  }
611 
613 
619  T& at(int index)
620  {
621  return const_cast<T&>(static_cast<const GridBase<T>* >(this)->at(index));
622  }
623 
625 
632  const T& at(const std::vector<double> &x) const
633  {
634  return at(GetIndices(x));
635  }
636 
638 
645  T& at(const std::vector<double> &x)
646  {
647  return at(GetIndices(x));
648  }
649 
651 
658  const T& at(double x) const
659  {
660  if (dimension_ != 1) {
661  throw std::invalid_argument("Only 1d-Grids can be accessed with a "
662  "single float as the specified point.");
663  }
664  return at({x});
665  }
666 
668 
674  T& at(double x)
675  {
676  return const_cast<T&>(static_cast<const GridBase<T>* >(this)->at(x));
677  }
678 
680 
684  const T& operator[](const std::vector<int> &indices) const
685  {
686  return at(indices);
687  }
688 
690 
694  T& operator[](const std::vector<int> &indices)
695  {
696  return at(indices);
697  }
698 
700 
710  template<typename R>
711  const T& operator[](std::initializer_list<R>&& x) const
712  {
713  return at(static_cast<std::vector<R> >(x));
714  }
715 
717 
727  template<typename R>
728  T& operator[](std::initializer_list<R>&& x)
729  {
730  return at(static_cast<std::vector<R> >(x));
731  }
732 
734 
740  const T& operator[](int index) const
741  {
742  return at(index);
743  }
744 
746 
752  T& operator[](int index)
753  {
754  return at(index);
755  }
756 
758 
762  const T& operator[](const std::vector<double> &x) const
763  {
764  return at(x);
765  }
766 
768 
772  T& operator[](const std::vector<double> &x)
773  {
774  return at(x);
775  }
776 
778 
784  const T& operator[](double x) const
785  {
786  return at(x);
787  }
788 
790 
796  T& operator[](double x)
797  {
798  return at(x);
799  }
800 
802  virtual ~GridBase() {}
803 };
804 
805 } // End namespace SSAGES
const std::vector< double > GetLower() const
Return the lower edges of the Grid.
Definition: GridBase.h:226
size_t size() const
Get the size of the internal storage vector.
Definition: GridBase.h:321
const T & at(const std::vector< int > &indices) const
Access Grid element read-only.
Definition: GridBase.h:541
const std::vector< bool > & GetPeriodic() const
Return the periodicity of the Grid.
Definition: GridBase.h:290
void syncGrid()
Sync the grid.
Definition: GridBase.h:141
std::vector< double > GetCoordinates(const std::vector< int > &indices)
! Return the distance to adjacent grid center points from given coordinates
Definition: GridBase.h:498
T & operator[](const std::vector< double > &x)
Access Grid element pertaining to a specific point per [] read-write.
Definition: GridBase.h:772
virtual ~GridBase()
Destructor.
Definition: GridBase.h:802
int GetNumPoints(size_t dim) const
Get the number of points for a specific dimension.
Definition: GridBase.h:212
T & at(int index)
Access 1d Grid by index, read-write.
Definition: GridBase.h:619
size_t GetDimension() const
Get the dimension.
Definition: GridBase.h:190
std::pair< std::vector< double >, std::vector< double > > edges_
Edges of the Grid in each dimension.
Definition: GridBase.h:51
STL namespace.
T & operator[](const std::vector< int > &indices)
Access Grid element per [] read-write.
Definition: GridBase.h:694
const T & operator[](double x) const
Access 1d-Grid via specific point, read-only.
Definition: GridBase.h:784
T & at(std::initializer_list< R > &&x)
Access Grid element via initializer list.
Definition: GridBase.h:591
T const * data() const
Get pointer to const of the internal data storage vector.
Definition: GridBase.h:345
const T & at(const std::vector< double > &x) const
Access Grid element pertaining to a specific point – read-only.
Definition: GridBase.h:632
T & at(const std::vector< int > &indices)
Access Grid element read/write.
Definition: GridBase.h:557
std::vector< T > data_
Internal storage of the data.
Definition: GridBase.h:42
std::vector< int > GetIndices(const std::vector< double > &x) const
Return the Grid indices for a given point.
Definition: GridBase.h:367
const T & at(int index) const
Access 1d Grid by index, read-only.
Definition: GridBase.h:603
const std::vector< double > GetUpper() const
Return the upper edges of the Grid.
Definition: GridBase.h:257
double GetLower(size_t dim) const
Get the lower edge for a specific dimension.
Definition: GridBase.h:242
bool GetPeriodic(size_t dim) const
Get the periodicity in a specific dimension.
Definition: GridBase.h:303
const std::vector< int > GetNumPoints() const
Get the number of points for all dimensions.
Definition: GridBase.h:200
std::vector< int > numPoints_
Number of points in each dimension.
Definition: GridBase.h:48
const T & operator[](const std::vector< double > &x) const
Access Grid element pertaining to a specific point per [] read-only.
Definition: GridBase.h:762
T & at(const std::vector< double > &x)
Access Grid element pertaining to a specific point – read/write.
Definition: GridBase.h:645
const T & operator[](std::initializer_list< R > &&x) const
Const access of Grid element via initializer list.
Definition: GridBase.h:711
std::vector< int > wrapIndices(const std::vector< int > &indices) const
Wrap the index around periodic boundaries.
Definition: GridBase.h:57
size_t dimension_
Dimension of the grid.
Definition: GridBase.h:45
virtual size_t mapTo1d(const std::vector< int > &indices) const =0
This function needs to be implemented by child classes.
double GetCoordinate(int index)
Return center point of 1d-grid.
Definition: GridBase.h:523
const T & operator[](int index) const
Access 1d-Grid per [] operator, read-only.
Definition: GridBase.h:740
T & operator[](int index)
Access 1d-Grid per [] operator, read-write.
Definition: GridBase.h:752
T & operator[](double x)
Access 1d-Grid via specific point, read-write.
Definition: GridBase.h:796
const T & at(std::initializer_list< R > &&x) const
Const access of Grid element via initializer list.
Definition: GridBase.h:574
T & at(double x)
Access 1d-Grid by point - read-write.
Definition: GridBase.h:674
double GetUpper(size_t dim) const
Get the upper edge for a specific dimension.
Definition: GridBase.h:274
const T & operator[](const std::vector< int > &indices) const
Access Grid element per [] read-only.
Definition: GridBase.h:684
T * data()
Get pointer to the internal data storage vector.
Definition: GridBase.h:333
GridBase(std::vector< int > numPoints, std::vector< double > lower, std::vector< double > upper, std::vector< bool > isPeriodic)
Constructor.
Definition: GridBase.h:99
T & operator[](std::initializer_list< R > &&x)
Access Grid element via initializer list.
Definition: GridBase.h:728
int GetIndex(double x) const
Return linear interpolation on a coordinate.
Definition: GridBase.h:467
std::vector< bool > isPeriodic_
Periodicity of the Grid.
Definition: GridBase.h:54
const T & at(double x) const
Access 1d-Grid by point - read-only.
Definition: GridBase.h:658
Base class for Grids.
Definition: GridBase.h:38