npstat is hosted by Hepforge, IPPP Durham
NPStat  5.10.0
mirrorWeight.hh
1 #ifndef NPSTAT_MIRRORWEIGHT_HH_
2 #define NPSTAT_MIRRORWEIGHT_HH_
3 
4 //======================================================================
5 // mirrorWeight.hh
6 //
7 // This is an internal header which is subject to change without notice.
8 // Application code should never call the functions declared/defined in
9 // this header directly.
10 //
11 // Author: I. Volobouev
12 //
13 // December 2011
14 //======================================================================
15 
16 #include <cassert>
17 #include <climits>
18 
19 #include "npstat/nm/ArrayND.hh"
20 
21 namespace npstat {
22  namespace Private {
23  template <typename T, unsigned StackLen, unsigned StackDim>
24  void mirrorWeightLoop(unsigned level, unsigned long idxFrom,
25  unsigned long idxTo, unsigned mirrorDirection,
26  const ArrayND<T,StackLen,StackDim>& from,
27  ArrayND<double,StackLen,StackDim>& to)
28  {
29  const bool positive(mirrorDirection & (1UL << level));
30  const unsigned imax = from.span(level);
31  if (level == from.rank() - 1)
32  {
33  const T* dfrom = from.data() + idxFrom;
34  double* dto = const_cast<double *>(to.data() + (idxTo+imax-1));
35  if (positive)
36  for (unsigned i=0; i<imax; ++i)
37  *dto++ = *dfrom++;
38  else
39  for (unsigned i=0; i<imax; ++i)
40  *dto-- = *dfrom++;
41  }
42  else
43  {
44  const unsigned long fromstride = from.strides()[level];
45  const unsigned long tostride = to.strides()[level];
46  idxTo += (imax - 1)*tostride;
47  for (unsigned i=0; i<imax; ++i)
48  {
49  mirrorWeightLoop(level+1, idxFrom, idxTo,
50  mirrorDirection, from, to);
51  idxFrom += fromstride;
52  if (positive)
53  idxTo += tostride;
54  else
55  idxTo -= tostride;
56  }
57  }
58  }
59 
60  //
61  // The "mirrorWeight" utility function turns the weight function scan in
62  // one hyperoctant into the complete scan over the whole support region.
63  // This can not be done by the ArrayND "multiMirror" method because the
64  // input hyperoctant includes the central grid points so that the result
65  // has an odd number of scanned points in each dimension.
66  //
67  template <typename T, unsigned StackLen, unsigned StackDim>
68  void mirrorWeight(const ArrayND<T,StackLen,StackDim>& from,
69  ArrayND<double,StackLen,StackDim>* to)
70  {
71  const unsigned dim = from.rank();
72  assert(dim);
73 #ifndef NDEBUG
74  const unsigned maxdim = CHAR_BIT*sizeof(unsigned long);
75  assert(dim < maxdim);
76 #endif
77  assert(to);
78  assert(dim == to->rank());
79  for (unsigned i=0; i<dim; ++i)
80  assert(to->span(i) == 2U*from.span(i) - 1U);
81  const unsigned long maxcycle = 1UL << dim;
82  for (unsigned long icycle=0UL; icycle<maxcycle; ++icycle)
83  mirrorWeightLoop(0U, 0UL, 0UL, icycle, from, *to);
84  }
85  }
86 }
87 
88 #endif // NPSTAT_MIRRORWEIGHT_HH_
Arbitrary-dimensional array template.
Definition: AbsArrayProjector.hh:14