QwGeant4
QweakSimFieldMap.hh
Go to the documentation of this file.
1 #ifndef QweakSimFieldMap_h
2 #define QweakSimFieldMap_h
3 
4 // System headers
5 #include <vector>
6 #include <fstream>
7 #include <iostream>
8 #include <stdexcept>
9 #include <stdint.h>
10 
11 // Geant4 headers
12 #include "G4ios.hh"
13 
14 /// Allowed interpolation methods
19 };
20 
21 // Define of the log stream and end line
22 #define mycout G4cout
23 #define mycerr G4cerr
24 #define myendl G4endl
25 
26 // Type of grid coordinates
27 typedef double coord_t;
28 
29 /**
30  * \class QweakSimFieldMap
31  * \ingroup QwAnalysis
32  * \brief A multi-dimensional grid of values with interpolation methods
33  *
34  * This class provides various interpolation methods on a multi-dimensional grid
35  * of multi-dimensional values. Linear interpolation and nearest-neighbor are
36  * implemented for all dimensions.
37  *
38  * The template arguments are the internal storage data type (defaults to float)
39  * and the number of dimensions of the stored data at each grid point (defaults
40  * to one). The dimension of the grid itself is set through the constructor.
41  * To describe a double vector field with 5 components on a 3-dimensional grid,
42  * you would write
43  * \code
44  * QwInterpolate<float,5> *field = new QweakSimFieldMap<float,5>(3);
45  * \endcode
46  *
47  * The minimum, maximum, and step size of the grid have to be known before
48  * the values are filled.
49  */
50 template <class value_t = float, unsigned int value_n = 1>
52 
53  public: // constructors and destructor
54 
55  /// Constructor with number of dimensions
56  QweakSimFieldMap(const unsigned int ndim = 1) {
57  SetDimensions(ndim);
59  };
60  /// Constructor with minimum, maximum, and step size
61  QweakSimFieldMap(const std::vector<coord_t>& min,
62  const std::vector<coord_t>& max,
63  const std::vector<coord_t>& step) {
64  SetDimensions(min.size());
66  SetMinimumMaximumStep(min,max,step);
67  };
68  /// Constructor with file name
69  QweakSimFieldMap(const std::string& filename) {
70  ReadBinaryFile(filename);
72  };
73  /// Destructor
74  virtual ~QweakSimFieldMap() { };
75 
76 
77  /// Debug level
79 
80  /// Set debug level
81  void SetDebugLevel(debug_t debug) { fDebug = debug; };
82 
83  private: // private member fields
84 
85  /// Number of dimensions in coordinates
86  unsigned int fNDim;
87 
88  /// Minimum in each dimension
89  std::vector<coord_t> fMin;
90  /// Maximum in each dimension
91  std::vector<coord_t> fMax;
92  /// Step size in each dimension
93  std::vector<coord_t> fStep;
94  /// Number of points in each dimension
95  std::vector<size_t> fSize;
96  /// Wrap around this coordinate
97  std::vector<size_t> fWrap;
98  /// Data reduction factor
99  std::vector<size_t> fRedux;
100 
101  /// Linear extent between neighbor points in each dimension (e.g. for the
102  /// least significant index this will be 1, for the next index the number
103  /// of points in the first index, etc...)
104  std::vector<size_t> fExtent;
105 
106  /// Table with pointers to arrays of values
107  std::vector<value_t> fValues[value_n];
108 
109  /// Number of values read in
110  unsigned int fCurrentEntries;
111  /// Maximum number of values
112  unsigned int fMaximumEntries;
113 
114  /// Interpolation method
116 
117  public: // public methods
118 
119  /// Set the number of coordinate dimensions and resize vectors
120  void SetDimensions(const unsigned int ndim) {
121  if (ndim == 0) {
122  mycerr << "QweakSimFieldMap::SetDimensions: "
123  << "Dimensions of the grid should be strictly positive!" << myendl;
124  return;
125  }
126  fNDim = ndim;
127  fMin.resize(fNDim); fMax.resize(fNDim); fStep.resize(fNDim); fWrap.resize(fNDim);
128  fSize.resize(fNDim); fExtent.resize(fNDim+1);
129  };
130  /// Set minimum, maximum, and step size to single values
131  void SetMinimumMaximumStep(const coord_t min, const coord_t max, const coord_t step) {
132  SetMinimumMaximumStep(std::vector<coord_t>(fNDim,min),
133  std::vector<coord_t>(fNDim,max),
134  std::vector<coord_t>(fNDim,step));
135  };
136  /// Set minimum, maximum, and step size to different values
137  void SetMinimumMaximumStep(const std::vector<coord_t>& min,
138  const std::vector<coord_t>& max,
139  const std::vector<coord_t>& step) {
140  // Set the dimensionality
141  if (min.size() != fNDim) SetDimensions(min.size());
142  // Check the dimensionality and assign boundaries and step size vectors
143  if (min.size() != fNDim || min.size() != fNDim || step.size() != fNDim) return;
144  fMin = min; fMax = max; fStep = step;
145  fExtent[0] = 1;
146  for (size_t i = 0; i < fMin.size(); i++) {
147  coord_t int_part; // safer to use modf than a direct static cast
148  coord_t frac_part = modf((fMax[i] - fMin[i]) / fStep[i] + 1.0, &int_part);
149  fSize[i] = static_cast<unsigned int>(int_part) + (frac_part > 0.5 ? 1 : 0);
150  fExtent[i+1] = fExtent[i] * fSize[i];
151  }
152  // Try resizing to allocate memory and initialize with zeroes
153  fMaximumEntries = fExtent[fNDim]; fCurrentEntries = 0; // no entries read yet
154  for (unsigned int i = 0; i < value_n; i++) {
155  try {
156  fValues[i].resize(fMaximumEntries,0);
157  } catch (std::bad_alloc) {
158  mycerr << "QweakSimFieldMap:SetMinimumMaximumStep: "
159  << "could not allocate memory for values!" << myendl;
160  }
161  }
162  };
163  /// Get minimum in dimension
164  coord_t GetMinimum(const unsigned int dim) const { return fMin.at(dim); };
165  /// Get maximum in dimension
166  coord_t GetMaximum(const unsigned int dim) const { return fMax.at(dim); };
167  /// Get minimum in dimension
168  coord_t GetStepSize(const unsigned int dim) const { return fStep.at(dim); };
169  /// Get the maximum number of entries
170  unsigned int GetMaximumEntries() const { return fMaximumEntries; };
171  /// Get the current number of entries
172  unsigned int GetCurrentEntries() const { return fCurrentEntries; };
173 
174  /// Get wrapping coordinate
175  unsigned int GetWrapCoordinate(const unsigned int dim) const
176  { return fWrap.at(dim); }
177  /// Set wrapping coordinate
178  void SetWrapCoordinate(const unsigned int dim, const size_t wrap = 1)
179  { fWrap.at(dim) = wrap; }
180  void SetWrapCoordinate(const std::vector<size_t>& wrap) {
181  if (wrap.size() != fNDim) return;
182  fWrap = wrap;
183  }
184 
185  /// Get data reduction factor
186  int GetDataReductionFactor(const unsigned int dim) const
187  { return fRedux.at(dim); }
188  /// Set data reduction factor
189  void SetDataReductionFactor(const unsigned int dim, const unsigned int redux)
190  { fRedux.at(dim) = redux; }
191  void SetDataReductionFactor(const unsigned int redux) {
192  for (unsigned int dim = 0; dim < fNDim; dim++)
193  SetDataReductionFactor(dim,redux);
194  }
195  void SetDataReductionFactor(const std::vector<unsigned int>& redux) {
196  if (redux.size() != fNDim) return;
197  fRedux = redux;
198  }
199 
200  /// Set the interpolation method
202  { fInterpolationMethod = method; };
203  /// Get the interpolation method
205  { return fInterpolationMethod; };
206 
207 
208  /// Print coverage map for all bins in one dimension
209  void PrintCoverage(const unsigned int dim) {
210  // Initialize coverage
211  unsigned int* cover = new unsigned int[fSize[dim]];
212  unsigned int* total = new unsigned int[fSize[dim]];
213  for (unsigned int i = 0; i < fSize[dim]; i++) {
214  cover[i] = 0;
215  total[i] = 0;
216  }
217  // Loop over all entries
218  value_t value[value_n];
219  unsigned int* cell_index = new unsigned int[fNDim];
220  for (unsigned int linear_index = 0; linear_index < fMaximumEntries; linear_index++) {
221  Cell(linear_index,cell_index);
222  total[cell_index[dim]]++;
223  if (Get(linear_index,value) && value[0] != 0)
224  cover[cell_index[dim]]++; // non-zero field
225  }
226  // Print coverage
227  coord_t* coord = new coord_t[fNDim];
228  for (size_t i = 0; i < fSize[dim]; i++) {
229  cell_index[dim] = i;
230  Coord(cell_index,coord);
231  mycout << "bin " << i << ", coord " << coord[dim] << ": "
232  << double(cover[i]) / double(total[i]) * 100 << "%"<< myendl;
233  }
234  delete[] cell_index;
235  delete[] coord;
236  delete[] cover;
237  delete[] total;
238  }
239 
240  /// Return true if the coordinate is in bounds
241  bool InBounds(const coord_t* coord) const {
242  return Check(coord);
243  };
244 
245 
246  /// \name Functions to write grid values
247  // @{
248  /// Set a single value at a coordinate (false if not possible)
249  bool Set(const coord_t& coord, const value_t& value) {
250  if (value_n != 1) {
251  mycerr << "QweakSimFieldMap::Set(coord_t&,value_t&): "
252  << "only valid for one-dimensional values." << myendl;
253  return false;
254  }
255  if (fNDim != 1) {
256  mycerr << "QweakSimFieldMap::Set(coord_t&,value_t&): "
257  << "only valid for one-dimensional grids." << myendl;
258  return false;
259  }
260  return Set(&coord, &value);
261  };
262  /// Set a single value at a coordinate (false if not possible)
263  bool Set(const coord_t* coord, const value_t& value) {
264  if (value_n != 1) {
265  mycerr << "QweakSimFieldMap::Set(coord_t*,value_t&): "
266  << "only valid for one-dimensional values." << myendl;
267  return false;
268  }
269  return Set(coord, &value);
270  };
271  /// Set a set of values at a coordinate (false if not possible)
272  bool Set(const coord_t* coord, const value_t* value) {
273  unsigned int* cell_index = new unsigned int[fNDim];
274  Nearest(coord, cell_index); // nearest cell
275  if (! Check(cell_index)) return false; // out of bounds
276  bool status = true;
277  bool written = false;
278  unsigned int linear_index;
279  for (unsigned int dim = 0; dim < fNDim; dim++) {
280  // skip dimensions that are not wrapped around
281  if (fWrap[dim] == 0) continue;
282  // FIXME only one wrapping coordinate correctly supported in Set()
283  if ((cell_index[dim] < fWrap[dim]) ||
284  (fSize[dim] - cell_index[dim] - 1 < fWrap[dim])) {
285  // there are equivalent grid points
286  for (size_t wrap = 0; wrap < fWrap[dim]; wrap++) {
287  // at the minimum
288  cell_index[dim] = wrap;
289  linear_index = Index(cell_index);
290  status &= Set(linear_index, value);
291  // at the maximum
292  cell_index[dim] = fSize[dim] - wrap - 1;
293  linear_index = Index(cell_index);
294  status &= Set(linear_index, value);
295  // set flag
296  written = true;
297  }
298  }
299  }
300  if (not written) {
301  // this is an unambiguous grid point
302  linear_index = Index(cell_index);
303  status &= Set(linear_index, value);
304  }
305  delete[] cell_index;
306  return status;
307  };
308 
309  /// Set a single value at a linearized index (false if not possible)
310  bool Set(const unsigned int linear_index, const value_t& value) {
311  if (value_n != 1) {
312  mycerr << "QweakSimFieldMap::Set(unsigned int,value_t&): "
313  << "only valid for one-dimensional values." << myendl;
314  return false;
315  }
316  return Set(linear_index, &value);
317  };
318  /// Set a set of values at a linearized index (false if not possible)
319  bool Set(const unsigned int linear_index, const value_t* value) {
320  if (! Check(linear_index)) return false; // out of bounds
321  for (unsigned int i = 0; i < value_n; i++)
322  fValues[i][linear_index] = value[i];
323  fCurrentEntries++;
324  return true;
325  };
326 
327  /// Set a single value at a grid point (false if out of bounds)
328  bool Set(const unsigned int* cell_index, const value_t& value) {
329  if (value_n != 1) {
330  mycerr << "QweakSimFieldMap::Set(unsigned int*,value_t&): "
331  << "only valid for one-dimensional values." << myendl;
332  return false;
333  }
334  return Set(cell_index, &value);
335  };
336  /// Set a set of values at a grid point (false if out of bounds)
337  bool Set(const unsigned int* cell_index, const value_t* value) {
338  if (! Check(cell_index)) return false; // out of bounds
339  return Set(Index(cell_index),value);
340  };
341  // @}
342 
343 
344  /// \name Functions to retrieve interpolated values
345  // @{
346  /// Get the interpolated value at coordinate (zero if out of bounds)
347  value_t GetValue(const coord_t& coord) const {
348  if (value_n != 1) {
349  mycerr << "QweakSimFieldMap::GetValue(coord_t&): "
350  << "only valid for one-dimensional fields." << myendl;
351  return false;
352  }
353  if (fNDim != 1) {
354  mycerr << "QweakSimFieldMap::GetValue(coord_t&): "
355  << "only valid for one-dimensional grids." << myendl;
356  return false;
357  }
358  value_t value;
359  if (GetValue(&coord, &value)) return value;
360  else return 0; // interpolation failed
361  };
362  /// Get the interpolated value at coordinate (zero if out of bounds)
363  value_t GetValue(const coord_t* coord) const {
364  if (value_n != 1) {
365  mycerr << "QweakSimFieldMap::GetValue(coord_t*): "
366  << "only valid for one-dimensional fields." << myendl;
367  return false;
368  }
369  value_t value;
370  if (GetValue(coord, &value)) return value;
371  else return 0; // interpolation failed
372  };
373  /// Get the interpolated value at coordinate (zero if out of bounds)
374  bool GetValue(const coord_t* coord, double& value) const {
375  if (value_n != 1) {
376  mycerr << "QweakSimFieldMap::GetValue(coord_t*,value_t&): "
377  << "only valid for one-dimensional fields." << myendl;
378  return false;
379  }
380  return GetValue(coord, &value);
381  };
382  /// Get the interpolated value at coordinate (zero if out of bounds)
383  bool GetValue(const coord_t* coord, double* value) const {
384  for (unsigned int i = 0; i < value_n; i++) value[i] = 0.0; // zero
385  if (! Check(coord)) {
386  if (fDebug >= kDebug) {
387  mycerr << "QweakSimFieldMap::GetValue(coord_t*,value_t*): "
388  << "coordinate is out of boundary." << myendl;
389  }
390  return false; // out of bounds
391  }
392  value_t result[value_n]; // we need a local copy of type value_t
393  switch (fInterpolationMethod) {
394  // return interpolation status
395  case kMultiLinear:
396  if (! Linear(coord, result)) {
397  if (fDebug >= kDebug) {
398  mycerr << "QweakSimFieldMap::GetValue(coord_t*,value_t*): "
399  << "linear interpolation failed." << myendl;
400  }
401  return false;
402  }
403  break;
404  case kNearestNeighbor:
405  if (! NearestNeighbor(coord, result)) {
406  if (fDebug >= kDebug) {
407  mycerr << "QweakSimFieldMap::GetValue(coord_t*,value_t*): "
408  << "nearest-neighbor interpolation failed." << myendl;
409  }
410  return false;
411  }
412  break;
413  default:
414  mycerr << "QweakSimFieldMap::GetValue(coord_t*,value_t*): "
415  << "unknown interpolation method." << myendl;
416  return false;
417  }
418  for (unsigned int i = 0; i < value_n; i++) value[i] = result[i]; // cast
419  return true;
420  };
421  // @}
422 
423 
424  /// \name File reading and writing
425  // @{
426  /// \brief Write the grid as text
427  bool WriteText(std::ostream& stream) const;
428  /// Write the grid to text file
429  bool WriteTextFile(const std::string& filename) const {
430  std::ofstream file(filename.c_str());
431  if (! file.is_open()) {
432  mycerr << "QweakSimFieldMap::WriteTextFile(string&): "
433  << "unable to open file " << filename << "." << myendl;
434  return false;
435  }
436  WriteText(file);
437  file.close();
438  return true;
439  };
440  /// Write the grid to screen
441  bool WriteTextScreen() const {
442  WriteText(mycout);
443  return true;
444  };
445  /// \brief Read the grid from text
446  bool ReadText(std::istream& stream);
447  /// Read the grid from text file
448  bool ReadTextFile(const std::string& filename) {
449  std::ifstream file(filename.c_str());
450  if (! file.is_open()) {
451  mycerr << "QweakSimFieldMap::ReadTextFile(string&): "
452  << "unable to open file " << filename << "." << myendl;
453  return false;
454  }
455  if (! ReadText(file)) {
456  mycerr << "QweakSimFieldMap::ReadTextFile(string&): "
457  << "unable to read text stream." << myendl;
458  return false;
459  }
460  file.close();
461  return true;
462  };
463  /// \brief Write the grid values to binary file
464  bool WriteBinaryFile(const std::string& filename) const;
465  /// \brief Read the grid values from binary file
466  bool ReadBinaryFile(const std::string& filename);
467  // @}
468 
469 
470  /// \name Indexing functions (publicly available and unchecked)
471  // @{
472  /// Return the linearized index based on the point coordinates (unchecked)
473  unsigned int Index(const coord_t* coord) const {
474  unsigned int* cell_index = new unsigned int[fNDim];
475  Cell(coord, cell_index);
476  unsigned int index = Index(cell_index);
477  delete[] cell_index;
478  return index;
479  };
480 
481  /// \brief Return the linearized index based on the cell indices (unchecked)
482  unsigned int Index(const unsigned int* cell_index) const;
483  /// \brief Return the linearized index based on the cell indices and offset (unchecked)
484  unsigned int Index(const unsigned int* cell_index, const unsigned int offset) const;
485 
486  /// \brief Return the cell index and local coordinates in one dimension (unchecked)
487  void Cell(const coord_t coord, unsigned int& cell_index, double& cell_local, const unsigned int dim) const;
488  /// \brief Return the cell index and local coordinates (unchecked)
489  void Cell(const coord_t* coord, unsigned int* cell_index, double* cell_local) const;
490  /// \brief Return the cell index (unchecked)
491  void Cell(const coord_t* coord, unsigned int* cell_index) const;
492  /// \brief Return the cell index based on the linearized index (unchecked)
493  void Cell(unsigned int linear_index, unsigned int* cell_index) const;
494 
495  /// \brief Return the coordinates based on the cell index (unchecked)
496  void Coord(const unsigned int* cell_index, coord_t* coord) const;
497  /// \brief Return the coordinates based on the linearized index (unchecked)
498  void Coord(const unsigned int linear_index, coord_t* coord) const;
499  // @}
500 
501  private: // private methods
502 
503  /// Return the cell index closest to the coordinate (could be above) (unchecked)
504  void Nearest(const coord_t* coord, unsigned int* cell_index) const {
505  double* cell_local = new double[fNDim];
506  Cell(coord, cell_index, cell_local);
507  // Loop over all dimensions and add one if larger than 0.5
508  for (unsigned int dim = 0; dim < fNDim; dim++)
509  if (cell_local[dim] > 0.5) cell_index[dim]++;
510  delete[] cell_local;
511  };
512 
513  /// \brief Linear interpolation (unchecked)
514  bool Linear(const coord_t* coord, value_t* value) const;
515  /// \brief Nearest-neighbor (unchecked)
516  bool NearestNeighbor(const coord_t* coord, value_t* value) const;
517 
518  /// \brief Check for boundaries with coordinate argument
519  bool Check(const coord_t* coord) const;
520  /// \brief Check for boundaries with cell index argument
521  bool Check(const unsigned int* cell_index) const;
522  /// \brief Check for boundaries with linearized index argument
523  bool Check(const unsigned int linear_index) const;
524 
525  /// Get a single value by cell index (unchecked)
526  value_t Get(const unsigned int* cell_index) const {
527  if (value_n != 1) {
528  mycerr << "QweakSimFieldMap::Get(unsigned int*): "
529  << "only valid for one-dimensional fields." << myendl;
530  return 0;
531  }
532  return fValues[0][Index(cell_index)];
533  };
534 
535  /// Get a single value by linearized index (unchecked)
536  value_t Get(const unsigned int index) const {
537  if (value_n != 1) {
538  mycerr << "QweakSimFieldMap::Get(unsigned int): "
539  << "only valid for one-dimensional fields." << myendl;
540  return 0;
541  }
542  return fValues[0][index];
543  };
544  /// Get a vector value by linearized index (unchecked)
545  bool Get(const unsigned int index, value_t* value) const {
546  for (unsigned int i = 0; i < value_n; i++)
547  value[i] = fValues[i][index];
548  return true;
549  };
550 
551 }; // class QweakSimFieldMap
552 
553 
554 
555 /**
556  * Perform the multi-dimensional linear interpolation (unchecked)
557  * @param coord Point coordinates
558  * @param value Interpolated value (reference)
559  */
560 template <class value_t, unsigned int value_n>
562  const coord_t* coord,
563  value_t* value) const
564 {
565  // Get cell and local normalized coordinates
566  unsigned int* cell_index = new unsigned int[fNDim];
567  double* cell_local = new double[fNDim];
568  Cell(coord, cell_index, cell_local);
569  // Initialize to zero
570  for (unsigned int i = 0; i < value_n; i++) value[i] = 0.0;
571  // Calculate the interpolated value
572  // by summing the 2^fNDim = 1 << fNDim neighbor points (1U is unsigned one)
573  for (unsigned int offset = 0; offset < (1U << fNDim); offset++) {
574  value_t neighbor[value_n];
575  if (! Get(Index(cell_index,offset), neighbor)) {
576  if (fDebug >= kDebug) {
577  mycout << "QweakSimFieldMap::Linear(coord_t*,value_t*): "
578  << "neighbor " << offset << " could not be found " << myendl;
579  }
580  return false;
581  }
582  // ... with appropriate weighting factors (1 - cell_local or cell_local)
583  double fac = 1.0;
584  for (unsigned int dim = 0; dim < fNDim; dim++)
585  fac *= ((offset >> dim) & 0x1)? cell_local[dim] : (1 - cell_local[dim]);
586  // ... for all vector components
587  for (unsigned int i = 0; i < value_n; i++)
588  value[i] += fac * neighbor[i];
589  }
590  delete[] cell_index;
591  delete[] cell_local;
592  return true;
593 }
594 
595 /**
596  * Perform the nearest-neighbor interpolation (unchecked)
597  * @param coord Point coordinates
598  * @param value Interpolated value (reference)
599  */
600 template <class value_t, unsigned int value_n>
602  const coord_t* coord,
603  value_t* value) const
604 {
605  // Get nearest cell
606  unsigned int* cell_index = new unsigned int[fNDim];
607  Nearest(coord, cell_index);
608  bool status = Get(Index(cell_index), value);
609  delete[] cell_index;
610  return status;
611 }
612 
613 /**
614  * Check whether the point is inside the valid region
615  * @param coord Point coordinates
616  * @return True if the coordinates are in the valid region
617  */
618 template <class value_t, unsigned int value_n>
619 inline bool QweakSimFieldMap<value_t,value_n>::Check(const coord_t* coord) const
620 {
621  for (unsigned int dim = 0; dim < fNDim; dim++) {
622  if (fWrap[dim] == 0 && (coord[dim] < fMin[dim] || coord[dim] > fMax[dim])) {
623  if (fDebug >= kDebug) {
624  mycout << "QweakSimFieldMap::Check(coord_t*): "
625  << "coord[" << dim << "] out of bounds "
626  << fMin[dim] << "-" << fMax[dim] << "." << myendl;
627  }
628  return false;
629  }
630  }
631  // Otherwise
632  return true;
633 }
634 
635 /**
636  * Check whether the cell index is inside the valid region
637  * @param cell_index Cell index
638  * @return True if the cell index is in the valid region
639  */
640 template <class value_t, unsigned int value_n>
641 inline bool QweakSimFieldMap<value_t,value_n>::Check(const unsigned int* cell_index) const
642 {
643  for (unsigned int dim = 0; dim < fNDim; dim++) {
644  if (cell_index[dim] >= fSize[dim]) {
645  if (fDebug >= kDebug) {
646  mycout << "QweakSimFieldMap::Check(unsigned int*): "
647  << "cell index[" << dim << "] " << cell_index[dim] << " larger "
648  << "than maximum " << fSize[dim] << "." << myendl;
649  }
650  return false;
651  }
652  }
653  // Otherwise
654  return true;
655 }
656 
657 /**
658  * Check whether the linearized index is inside the valid region
659  * @param linear_index Linearized index
660  * @return True if the cell index is in the valid region
661  */
662 template <class value_t, unsigned int value_n>
663 inline bool QweakSimFieldMap<value_t,value_n>::Check(const unsigned int linear_index) const
664 {
665  if (linear_index >= fMaximumEntries) {
666  if (fDebug >= kDebug) {
667  mycout << "QweakSimFieldMap::Check(unsigned int): "
668  << "linear index " << linear_index << " larger "
669  << "than maximum " << fMaximumEntries << "." << myendl;
670  }
671  return false;
672  }
673  // Otherwise
674  return true;
675 }
676 
677 /**
678  * Return the linearized index based on the cell indices (unchecked)
679  * @param cell_index Index in each dimension
680  * @return Linearized index
681  */
682 template <class value_t, unsigned int value_n>
684  const unsigned int* cell_index) const
685 {
686  unsigned int linear_index = 0;
687  // Loop over dimensions
688  for (unsigned int dim = 0; dim < fNDim; dim++)
689  // ... and use the stored extents for every dimensions
690  linear_index += fExtent[dim] * cell_index[dim];
691  return linear_index;
692 }
693 
694 /**
695  * Return the linearized index based on the cell indices and offset (unchecked)
696  * @param cell_index Index in each dimension
697  * @param pattern Bit pattern with offsets in each dimension
698  * @return Linearized index
699  */
700 template <class value_t, unsigned int value_n>
702  const unsigned int* cell_index,
703  const unsigned int pattern) const
704 {
705  unsigned int linear_index = 0;
706  // Loop over dimensions
707  for (unsigned int dim = 0; dim < fNDim; dim++) {
708  // Add offset of one based on the bit at position dim
709  unsigned int offset = (pattern >> dim) & 0x1;
710  linear_index += fExtent[dim] * (cell_index[dim] + offset);
711  }
712  return linear_index;
713 }
714 
715 /**
716  * Return the cell index and local coordinates in one dimension (unchecked)
717  * @param coord Point coordinate in one dimension
718  * @param cell_index Cell index of the point (reference)
719  * @param cell_local Local coordinates in cell (reference)
720  * @param dim Dimension
721  */
722 template <class value_t, unsigned int value_n>
724  const coord_t coord,
725  unsigned int &cell_index,
726  double &cell_local,
727  const unsigned int dim) const
728 {
729  // Normalize coordinate (positive, with integers on grid points)
730  double norm_coord = (coord - fMin[dim]) / fStep[dim];
731  // Split in fractional and integer part
732  double frac_part, int_part;
733  frac_part = modf(norm_coord, &int_part);
734  cell_local = frac_part;
735  cell_index = static_cast<int>(int_part); // cast to integer
736  // Wrap index
737  if (fWrap[dim] > 0) cell_index %= (fSize[dim] - fWrap[dim]);
738 }
739 
740 /**
741  * Return the cell index and local coordinates (unchecked)
742  * @param coord Point coordinates
743  * @param cell_index Cell index of the point (reference)
744  * @param cell_local Local coordinates in cell (reference)
745  */
746 template <class value_t, unsigned int value_n>
748  const coord_t* coord,
749  unsigned int* cell_index,
750  double* cell_local) const
751 {
752  // Get cell index and local coordinates in each dimension
753  for (unsigned int dim = 0; dim < fNDim; dim++)
754  Cell(coord[dim], cell_index[dim], cell_local[dim], dim);
755 }
756 
757 /**
758  * Return the cell index (unchecked)
759  * @param coord Point coordinates
760  * @param cell_index Cell index of the point (reference)
761  */
762 template <class value_t, unsigned int value_n>
764  const coord_t* coord,
765  unsigned int* cell_index) const
766 {
767  // Get cell index and ignore local coordinates
768  double* cell_local = new double[fNDim];
769  Cell(coord, cell_index, cell_local);
770  delete[] cell_local;
771 }
772 
773 /**
774  * Return the cell index based on the linearized index (unchecked)
775  * @param linear_index Linearized index
776  * @param cell_index Cell index (reference)
777  */
778 template <class value_t, unsigned int value_n>
780  unsigned int linear_index,
781  unsigned int* cell_index) const
782 {
783  // This does not work with unsigned int, because that is always >= 0 and wraps
784  for (int dim = fNDim-1; dim >= 0; dim--) {
785  cell_index[dim] = linear_index / fExtent[dim];
786  linear_index -= cell_index[dim] * fExtent[dim];
787  }
788 }
789 
790 /**
791  * Return the coordinates based on the cell index (unchecked)
792  * @param cell_index Cell index
793  * @param coord Point coordinates (reference)
794  */
795 template <class value_t, unsigned int value_n>
797  const unsigned int* cell_index,
798  coord_t* coord) const
799 {
800  for (unsigned int dim = 0; dim < fNDim; dim++)
801  coord[dim] = fMin[dim] + cell_index[dim] * fStep[dim];
802 }
803 
804 /**
805  * Return the coordinates based on the linearized index (unchecked)
806  * @param linear_index Linearized index
807  * @param coord Point coordinates (reference)
808  */
809 template <class value_t, unsigned int value_n>
811  const unsigned int linear_index,
812  coord_t* coord) const
813 {
814  unsigned int* cell_index = new unsigned int[fNDim];
815  Cell(linear_index,cell_index);
816  Coord(cell_index,coord);
817  delete[] cell_index;
818 }
819 
820 /**
821  * Write the grid values to a text stream
822  * @param stream Output stream
823  */
824 template <class value_t, unsigned int value_n>
825 inline bool QweakSimFieldMap<value_t,value_n>::WriteText(std::ostream& stream) const
826 {
827  // Informational message
828  mycout << "Writing text stream: ";
829  // Write the dimensions
830  stream << fNDim << "\t" << value_n << std::endl;
831  // Write the grid parameters
832  for (unsigned int dim = 0; dim < fNDim; dim++)
833  stream << dim << "\t" << fMin[dim] << "\t" << fMax[dim] << "\t" << fStep[dim] << std::endl;
834  // Write the values
835  stream << fValues[0].size() << std::endl;
836  unsigned int entries = fValues[0].size();
837  for (unsigned int index = 0; index < entries; index++) {
838  // Write values
839  for (unsigned int i = 0; i < value_n; i++) {
840  stream << fValues[i][index] << "\t";
841  }
842  stream << std::endl;
843  // Progress bar
844  if (index % (entries / 10) == 0)
845  mycout << index / (entries / 100) << "%" << std::flush;
846  if (index % (entries / 10) != 0
847  && index % (entries / 40) == 0)
848  mycout << "." << std::flush;
849  }
850  stream << "end" << std::endl;
851  mycout << myendl;
852  return true;
853 }
854 
855 /**
856  * Read the grid values from a text stream
857  * @param stream Input stream
858  * @return True if successfully read all values
859  */
860 template <class value_t, unsigned int value_n>
861 inline bool QweakSimFieldMap<value_t,value_n>::ReadText(std::istream& stream)
862 {
863  // Informational message
864  mycout << "Reading text stream: ";
865  // Read the dimensions
866  unsigned int n;
867  stream >> fNDim >> n;
868  if (n != value_n) {
869  mycerr << "QweakSimFieldMap::ReadText(istream&): "
870  << "incompatible field dimension in stream." << myendl;
871  return false;
872  }
873  SetDimensions(fNDim);
874  // Read the grid parameters
875  for (unsigned int dim = 0; dim < fNDim; dim++)
876  stream >> dim >> fMin[dim] >> fMax[dim] >> fStep[dim];
877  SetMinimumMaximumStep(fMin, fMax, fStep);
878  // Read the grid values
879  unsigned int entries;
880  stream >> entries;
881  for (unsigned int index = 0; index < entries; index++) {
882  // Read values
883  for (unsigned int i = 0; i < value_n; i++) {
884  stream >> fValues[i][index];
885  }
886  // Progress bar
887  if (index % (entries / 10) == 0)
888  mycout << index / (entries / 100) << "%" << std::flush;
889  if (index % (entries / 10) != 0
890  && index % (entries / 40) == 0)
891  mycout << "." << std::flush;
892  }
893  // Check for end of file
894  std::string end;
895  stream >> end;
896  // Informational message
897  mycout << myendl;
898  if (end == "end") return true;
899  else {
900  mycerr << "QweakSimFieldMap::ReadText(istream&): "
901  << "stream did not end with keyword 'end'." << myendl;
902  return false;
903  }
904 }
905 
906 /**
907  * Write the grid values to binary file (should be 64-bit safe, untested)
908  *
909  * Integer data types can be stored differently on 32-bit and 64-bit systems.
910  * Fixed-length types uint32_t and u_int32_t are provided in stdint.h and
911  * sys/types.h, respectively. The floating point types float and double will
912  * always have a length of 32 and 64 bit, per the IEEE convention.
913  *
914  * @param filename File name
915  * @return True if written successfully
916  */
917 template <class value_t, unsigned int value_n>
918 inline bool QweakSimFieldMap<value_t,value_n>::WriteBinaryFile(const std::string& filename) const
919 {
920  std::ofstream file(filename.c_str(), std::ios::binary);
921  if (! file.is_open()) {
922  mycerr << "QweakSimFieldMap::WriteBinaryFile(string&): "
923  << "unable to open file " << filename << "." << myendl;
924  return false;
925  }
926  // Informational message
927  mycout << "Writing binary file: ";
928  // Write template parameters
929  uint32_t n = value_n; // uint32_t has length of 32 bits on any system
930  uint32_t size = sizeof(value_t);
931  file.write(reinterpret_cast<const char*>(&n), sizeof(n));
932  file.write(reinterpret_cast<const char*>(&size), sizeof(size));
933  // Write grid parameters
934  uint32_t ndim = fNDim;
935  file.write(reinterpret_cast<const char*>(&ndim), sizeof(ndim));
936  for (unsigned int dim = 0; dim < fNDim; dim++) {
937  file.write(reinterpret_cast<const char*>(&fMin[dim]),sizeof(fMin[dim]));
938  file.write(reinterpret_cast<const char*>(&fMax[dim]),sizeof(fMax[dim]));
939  file.write(reinterpret_cast<const char*>(&fStep[dim]),sizeof(fStep[dim]));
940  }
941  uint32_t entries = fValues[0].size();
942  file.write(reinterpret_cast<const char*>(&entries),sizeof(entries));
943  for (unsigned int index = 0; index < entries; index++) {
944  // Write values
945  for (unsigned int i = 0; i < value_n; i++) {
946  value_t value = fValues[i][index];
947  file.write(reinterpret_cast<const char*>(&value),sizeof(value));
948  }
949  // Progress bar
950  if (index % (entries / 10) == 0)
951  mycout << index / (entries / 100) << "%" << std::flush;
952  if (index % (entries / 10) != 0
953  && index % (entries / 40) == 0)
954  mycout << "." << std::flush;
955  }
956  mycout << myendl;
957  file.close();
958  return true;
959 }
960 
961 /**
962  * Read the grid values from binary file (should be 64-bit safe, untested)
963  * @param filename File name
964  * @return True if read successfully
965  */
966 template <class value_t, unsigned int value_n>
967 inline bool QweakSimFieldMap<value_t,value_n>::ReadBinaryFile(const std::string& filename)
968 {
969  std::ifstream file(filename.c_str(), std::ios::binary);
970  if (! file.is_open()) {
971  mycerr << "QweakSimFieldMap::ReadBinaryFile(string&): "
972  << "unable to open file " << filename << "." << myendl;
973  return false;
974  }
975  // Informational message
976  mycout << "Reading binary file: ";
977  // Go to end and store length (could also use std::ios::ate)
978  file.seekg(0, std::ios::end);
979  // Go to begin and start reading template parameters
980  file.seekg(0, std::ios::beg);
981  uint32_t n, size;
982  file.read(reinterpret_cast<char*>(&n), sizeof(n));
983  file.read(reinterpret_cast<char*>(&size), sizeof(size));
984  if (n != value_n) {
985  mycerr << "QweakSimFieldMap::ReadBinaryFile(string&): "
986  << "incompatible field dimension in stream." << myendl;
987  return false;
988  }
989  if (size != sizeof(value_t)) {
990  mycerr << "QweakSimFieldMap::ReadBinaryFile(string&): "
991  << "incompatible precision of field values.." << myendl;
992  return false;
993  }
994  // Read grid parameters
995  uint32_t ndim;
996  file.read(reinterpret_cast<char*>(&ndim), sizeof(ndim));
997  SetDimensions(ndim);
998  // Read sizes
999  for (unsigned int dim = 0; dim < fNDim; dim++) {
1000  file.read(reinterpret_cast<char*>(&fMin[dim]),sizeof(fMin[dim]));
1001  file.read(reinterpret_cast<char*>(&fMax[dim]),sizeof(fMax[dim]));
1002  file.read(reinterpret_cast<char*>(&fStep[dim]),sizeof(fStep[dim]));
1003  }
1004  SetMinimumMaximumStep(fMin, fMax, fStep); // calculates fMaximumEntries
1005  // Read values
1006  uint32_t maximum_entries;
1007  file.read(reinterpret_cast<char*>(&maximum_entries),sizeof(maximum_entries));
1008  if (maximum_entries != fMaximumEntries) {
1009  mycerr << "QweakSimFieldMap::ReadBinaryFile(string&): "
1010  << "expected number of entries." << myendl;
1011  return false;
1012  }
1013  int value_size = sizeof(value_t);
1014  unsigned int entries = fValues[0].size();
1015  for (unsigned int index = 0; index < entries; index++) {
1016  // Read values
1017  for (unsigned int i = 0; i < value_n; i++) {
1018  file.read((char*)(&fValues[i][index]),value_size);
1019  }
1020  // Progress bar
1021  if (index % (entries / 10) == 0)
1022  mycout << index / (entries / 100) << "%" << std::flush;
1023  if (index % (entries / 10) != 0
1024  && index % (entries / 40) == 0)
1025  mycout << "." << std::flush;
1026  }
1027  mycout << myendl;
1028  file.close();
1029  return true;
1030 }
1031 
1032 
1033 #endif // QweakSimFieldMap_h
1034 
EInterpolationMethod GetInterpolationMethod() const
Get the interpolation method.
value_t GetValue(const coord_t &coord) const
Get the interpolated value at coordinate (zero if out of bounds)
int GetDataReductionFactor(const unsigned int dim) const
Get data reduction factor.
std::vector< coord_t > fMax
Maximum in each dimension.
bool Set(const unsigned int *cell_index, const value_t *value)
Set a set of values at a grid point (false if out of bounds)
virtual ~QweakSimFieldMap()
Destructor.
bool InBounds(const coord_t *coord) const
Return true if the coordinate is in bounds.
bool WriteTextScreen() const
Write the grid to screen.
bool NearestNeighbor(const coord_t *coord, value_t *value) const
Nearest-neighbor (unchecked)
value_t Get(const unsigned int index) const
Get a single value by linearized index (unchecked)
value_t Get(const unsigned int *cell_index) const
Get a single value by cell index (unchecked)
coord_t GetMinimum(const unsigned int dim) const
Get minimum in dimension.
enum QweakSimFieldMap::debug_t fDebug
void SetMinimumMaximumStep(const coord_t min, const coord_t max, const coord_t step)
Set minimum, maximum, and step size to single values.
bool ReadText(std::istream &stream)
Read the grid from text.
bool ReadTextFile(const std::string &filename)
Read the grid from text file.
std::vector< value_t > fValues[value_n]
Table with pointers to arrays of values.
unsigned int GetCurrentEntries() const
Get the current number of entries.
double coord_t
std::vector< size_t > fExtent
A multi-dimensional grid of values with interpolation methods.
void Coord(const unsigned int *cell_index, coord_t *coord) const
Return the coordinates based on the cell index (unchecked)
bool WriteBinaryFile(const std::string &filename) const
Write the grid values to binary file.
bool WriteTextFile(const std::string &filename) const
Write the grid to text file.
void SetDataReductionFactor(const unsigned int dim, const unsigned int redux)
Set data reduction factor.
unsigned int GetMaximumEntries() const
Get the maximum number of entries.
bool Check(const coord_t *coord) const
Check for boundaries with coordinate argument.
EInterpolationMethod fInterpolationMethod
Interpolation method.
bool Set(const coord_t *coord, const value_t *value)
Set a set of values at a coordinate (false if not possible)
unsigned int fMaximumEntries
Maximum number of values.
debug_t
Debug level.
bool Set(const unsigned int *cell_index, const value_t &value)
Set a single value at a grid point (false if out of bounds)
value_t GetValue(const coord_t *coord) const
Get the interpolated value at coordinate (zero if out of bounds)
unsigned int GetWrapCoordinate(const unsigned int dim) const
Get wrapping coordinate.
coord_t GetMaximum(const unsigned int dim) const
Get maximum in dimension.
bool Get(const unsigned int index, value_t *value) const
Get a vector value by linearized index (unchecked)
QweakSimFieldMap(const std::string &filename)
Constructor with file name.
std::vector< coord_t > fStep
Step size in each dimension.
std::vector< size_t > fRedux
Data reduction factor.
unsigned int fNDim
Number of dimensions in coordinates.
void Cell(const coord_t coord, unsigned int &cell_index, double &cell_local, const unsigned int dim) const
Return the cell index and local coordinates in one dimension (unchecked)
void PrintCoverage(const unsigned int dim)
Print coverage map for all bins in one dimension.
std::vector< size_t > fWrap
Wrap around this coordinate.
void SetMinimumMaximumStep(const std::vector< coord_t > &min, const std::vector< coord_t > &max, const std::vector< coord_t > &step)
Set minimum, maximum, and step size to different values.
bool Set(const coord_t &coord, const value_t &value)
Set a single value at a coordinate (false if not possible)
std::vector< size_t > fSize
Number of points in each dimension.
bool Linear(const coord_t *coord, value_t *value) const
Linear interpolation (unchecked)
#define mycerr
std::vector< coord_t > fMin
Minimum in each dimension.
#define mycout
coord_t GetStepSize(const unsigned int dim) const
Get minimum in dimension.
bool Set(const coord_t *coord, const value_t &value)
Set a single value at a coordinate (false if not possible)
bool Set(const unsigned int linear_index, const value_t &value)
Set a single value at a linearized index (false if not possible)
void SetWrapCoordinate(const unsigned int dim, const size_t wrap=1)
Set wrapping coordinate.
unsigned int fCurrentEntries
Number of values read in.
bool GetValue(const coord_t *coord, double *value) const
Get the interpolated value at coordinate (zero if out of bounds)
void SetWrapCoordinate(const std::vector< size_t > &wrap)
#define myendl
bool ReadBinaryFile(const std::string &filename)
Read the grid values from binary file.
void SetDimensions(const unsigned int ndim)
Set the number of coordinate dimensions and resize vectors.
void SetDataReductionFactor(const std::vector< unsigned int > &redux)
void SetDataReductionFactor(const unsigned int redux)
QweakSimFieldMap(const std::vector< coord_t > &min, const std::vector< coord_t > &max, const std::vector< coord_t > &step)
Constructor with minimum, maximum, and step size.
void SetDebugLevel(debug_t debug)
Set debug level.
bool GetValue(const coord_t *coord, double &value) const
Get the interpolated value at coordinate (zero if out of bounds)
EInterpolationMethod
Allowed interpolation methods.
unsigned int Index(const coord_t *coord) const
Return the linearized index based on the point coordinates (unchecked)
bool Set(const unsigned int linear_index, const value_t *value)
Set a set of values at a linearized index (false if not possible)
QweakSimFieldMap(const unsigned int ndim=1)
Constructor with number of dimensions.
void Nearest(const coord_t *coord, unsigned int *cell_index) const
Return the cell index closest to the coordinate (could be above) (unchecked)
void SetInterpolationMethod(const EInterpolationMethod method)
Set the interpolation method.
bool WriteText(std::ostream &stream) const
Write the grid as text.