npstat is hosted by Hepforge, IPPP Durham
NPStat  5.10.0
LocalMultiFilter1D.hh
Go to the documentation of this file.
1 #ifndef NPSTAT_LOCALMULTIFILTER1D_HH_
2 #define NPSTAT_LOCALMULTIFILTER1D_HH_
3 
4 /*!
5 // \file LocalMultiFilter1D.hh
6 //
7 // \brief Local polynomial filtering (regression) on 1-d equidistant grids
8 //
9 // Generates a collection of filters. Each filter in this collection
10 // corresponds to a polynomial of a certain degree from the orthogonal
11 // polynomial set, with degrees varying from 0 to the specified maximum.
12 // Intended mainly for numerical exploration of such filters.
13 //
14 // Author: I. Volobouev
15 //
16 // April 2015
17 */
18 
19 #include <cassert>
20 #include <stdexcept>
21 
22 #include "geners/CPP11_auto_ptr.hh"
23 #include "geners/ClassId.hh"
24 
25 #include "npstat/nm/Matrix.hh"
26 
29 
30 #ifdef SWIG
31 #include "npstat/wrap/arrayNDToNumpy.hh"
32 #include "npstat/wrap/numpyArrayUtils.hh"
33 #endif // SWIG
34 
35 namespace npstat {
37  {
38  public:
39  /**
40  // Main constructor. The arguments are as follows:
41  //
42  // maxDegree -- Maximum degree of the polynomials.
43  //
44  // filterBuilder -- An instance of a class which actually builds
45  // the filters when requested. This builder
46  // is used only inside the constructor.
47  //
48  // dataLen -- The length of the data arrays which will be
49  // used with this filter (this info is needed in
50  // order to take into account the boundary effects).
51  */
53  const OrthoPolyFilter1DBuilder& filterBuilder,
54  unsigned dataLen);
55 
57  LocalMultiFilter1D& operator=(const LocalMultiFilter1D&);
58 
60 
61  bool operator==(const LocalMultiFilter1D& r) const;
62  inline bool operator!=(const LocalMultiFilter1D& r) const
63  {return !(*this == r);}
64 
65  //@{
66  /** Inspect object properties */
67  inline unsigned maxDegree() const {return maxDegree_;}
68  inline unsigned dataLen() const {return nbins_;}
69  //@}
70 
71  /** Get the filter coefficients for the given bin */
72  const PolyFilter1D& getFilter(unsigned degree, unsigned binNumber) const;
73 
74  /**
75  // This method performs the filtering. "dataLen",
76  // which is the length of both "in" and "out" arrays,
77  // must be the same as the one in the constructor.
78  */
79  template <typename Tin, typename Tout>
80  void filter(unsigned degree, const Tin* in,
81  unsigned dataLen, Tout* out) const;
82 
83  /**
84  // A diffent filtering method in which the shapes of the
85  // kernels are determined by the positions of the "sources"
86  // (i.e., sample points) instead of the positions at which
87  // the density (or response) is estimated.
88  */
89  template <typename Tin, typename Tout>
90  void convolve(unsigned degree, const Tin* in,
91  unsigned dataLen, Tout* out) const;
92 
93  /**
94  // Generate the complete (non-sparse) representation of the filter.
95  // It will be a generalized stochastic matrix (each row sums to 1).
96  */
97  Matrix<double> getFilterMatrix(unsigned degree) const;
98 
99  // Methods needed for I/O
100  inline gs::ClassId classId() const {return gs::ClassId(*this);}
101  bool write(std::ostream& os) const;
102 
103  static inline const char* classname()
104  {return "npstat::LocalMultiFilter1D";}
105  static inline unsigned version() {return 1;}
106  static LocalMultiFilter1D* read(const gs::ClassId& id, std::istream& in);
107 
108  private:
110 
111  mutable std::vector<long double> sumBuffer_;
112  std::vector<PolyFilter1D*> unique_;
113  std::vector<PolyFilter1D**> bins_;
114  unsigned nbins_;
115  unsigned maxDegree_;
116 
117  template <typename T>
118  double convolute(const T* data, const PolyFilter1D* filter) const;
119 
120  void addWeightedFilter(long double w, unsigned degree,
121  unsigned binNum) const;
122 
123  void clearSumBuffer() const;
124 
125  void releaseMem();
126  void copyOtherData(const LocalMultiFilter1D&);
127 
128  // The following function will assume that is was
129  // given a matrix filled with zeros as an argument
130  void fillFilterMatrix(unsigned degree, Matrix<double>* fm) const;
131 
132 #ifdef SWIG
133  public:
134  inline PyObject* filter_2(unsigned degree,
135  const double* in, unsigned dataLen,
136  const bool makeNonNegative = false) const
137  {
138  if (dataLen != nbins_) throw std::invalid_argument(
139  "In npstat::LocalPolyFilter1D::filter_2: "
140  "incompatible data length");
141  assert(in);
142 
143  const int typenum = NumpyTypecode<double>::code;
144  npy_intp sh = dataLen;
145  PyObject* xarr = PyArray_SimpleNew(1, &sh, typenum);
146  if (xarr)
147  {
148  double* out = (double*)PyArray_DATA((PyArrayObject*)xarr);
149  try
150  {
151  this->filter(degree, in, dataLen, out);
152  }
153  catch (const std::exception& e)
154  {
155  Py_DECREF(xarr);
156  throw e;
157  }
158  if (makeNonNegative)
159  for (unsigned i=0; i<dataLen; ++i)
160  if (out[i] < 0.0)
161  out[i] = 0.0;
162  }
163  return xarr;
164  }
165 
166  inline PyObject* convolve_2(unsigned degree,
167  const double* in, unsigned dataLen,
168  const bool makeNonNegative = false) const
169  {
170  if (dataLen != nbins_) throw std::invalid_argument(
171  "In npstat::LocalPolyFilter1D::convolve_2: "
172  "incompatible data length");
173  assert(in);
174 
175  const int typenum = NumpyTypecode<double>::code;
176  npy_intp sh = dataLen;
177  PyObject* xarr = PyArray_SimpleNew(1, &sh, typenum);
178  if (xarr)
179  {
180  double* out = (double*)PyArray_DATA((PyArrayObject*)xarr);
181  try
182  {
183  this->convolve(degree, in, dataLen, out);
184  }
185  catch (const std::exception& e)
186  {
187  Py_DECREF(xarr);
188  throw e;
189  }
190  if (makeNonNegative)
191  for (unsigned i=0; i<dataLen; ++i)
192  if (out[i] < 0.0)
193  out[i] = 0.0;
194  }
195  return xarr;
196  }
197 
198  inline void filterIntoBuffer(unsigned degree,
199  const double* in, unsigned dataLen,
200  double* out, unsigned outLen,
201  const bool makeNonNegative = false) const
202  {
203  if (dataLen != outLen) throw std::invalid_argument(
204  "In npstat::LocalMultiFilter1D::filterIntoBuffer: "
205  "incompatible array sizes");
206  this->filter(degree, in, dataLen, out);
207  if (makeNonNegative)
208  for (unsigned i=0; i<dataLen; ++i)
209  if (out[i] < 0.0)
210  out[i] = 0.0;
211  }
212 
213  inline void convolveIntoBuffer(unsigned degree,
214  const double* in, unsigned dataLen,
215  double* out, unsigned outLen,
216  const bool makeNonNegative = false) const
217  {
218  if (dataLen != outLen) throw std::invalid_argument(
219  "In npstat::LocalMultiFilter1D::convolveIntoBuffer: "
220  "incompatible array sizes");
221  this->convolve(degree, in, dataLen, out);
222  if (makeNonNegative)
223  for (unsigned i=0; i<dataLen; ++i)
224  if (out[i] < 0.0)
225  out[i] = 0.0;
226  }
227 #endif // SWIG
228  };
229 
230  /**
231  // The utility function that generates the most common filters using
232  // kernels from the symmetric beta family. The arguments are as follows:
233  //
234  // m -- Choose the kernel from the symmetric beta family
235  // proportional to (1 - x^2)^m. If m is negative,
236  // Gaussian kernel will be used, truncated at +- 12 sigma.
237  //
238  // bandwidth -- Kernel bandwidth.
239  //
240  // maxDegree -- Maximum degree of the LOrPE polynomial.
241  //
242  // numberOfGridPoints -- Length of the data array to be used with this
243  // filter (typically, number of histogram bins
244  // and such)
245  //
246  // xmin, xmax -- Data grid limits (as in a histogram).
247  //
248  // boundaryMethod -- Method for handling the weight function at
249  // the boundary of the density support region.
250  //
251  // exclusionMask -- If provided, array with numberOfGridPoints elements.
252  // If an element of this array is not 0, corresponding
253  // data point will be excluded from the filtering process.
254  //
255  // excludeCentralPoint -- If "true", the weight will be set to 0 for
256  // the central bin of the filter. This can be useful in
257  // some cross validation scenarios.
258  */
259  CPP11_auto_ptr<LocalMultiFilter1D> symbetaMultiFilter1D(
260  int m, double bandwidth, unsigned maxDegree,
261  unsigned numberOfGridPoints, double xmin, double xmax,
262  const BoundaryHandling& boundaryMethod,
263  const unsigned char* exclusionMask = 0,
264  bool excludeCentralPoint = false);
265 
266 #ifdef SWIG
267  inline LocalMultiFilter1D* symbetaMultiFilter1D_2(
268  int m, double bandwidth, unsigned maxDegree,
269  unsigned numberOfGridPoints, double xmin, double xmax,
270  const BoundaryHandling& boundaryMethod,
271  const unsigned char* exclusionMask, unsigned maskLength,
272  bool excludeCentralPoint = false)
273  {
274  const unsigned char* mask = 0;
275  if (maskLength)
276  {
277  assert(exclusionMask);
278  mask = exclusionMask;
279  if (maskLength != numberOfGridPoints) throw std::invalid_argument(
280  "In npstat::symbetaMultiFilter1D_2: incompatible mask length");
281  }
282  CPP11_auto_ptr<LocalMultiFilter1D> ptr = symbetaMultiFilter1D(
283  m, bandwidth, maxDegree, numberOfGridPoints, xmin, xmax,
284  boundaryMethod, mask, excludeCentralPoint);
285  return ptr.release();
286  }
287 #endif // SWIG
288 }
289 
290 #include "npstat/stat/LocalMultiFilter1D.icc"
291 
292 #endif // NPSTAT_LOCALMULTIFILTER1D_HH_
Abstract interface for building local polynomial filter weights in 1-d.
API for LOrPE boundary handling methods.
Template matrix class.
Definition: BoundaryHandling.hh:21
Definition: LocalMultiFilter1D.hh:37
void filter(unsigned degree, const Tin *in, unsigned dataLen, Tout *out) const
unsigned maxDegree() const
Definition: LocalMultiFilter1D.hh:67
void convolve(unsigned degree, const Tin *in, unsigned dataLen, Tout *out) const
LocalMultiFilter1D(unsigned maxDegree, const OrthoPolyFilter1DBuilder &filterBuilder, unsigned dataLen)
const PolyFilter1D & getFilter(unsigned degree, unsigned binNumber) const
Matrix< double > getFilterMatrix(unsigned degree) const
Definition: AbsFilter1DBuilder.hh:34
Definition: AbsArrayProjector.hh:14
CPP11_auto_ptr< LocalMultiFilter1D > symbetaMultiFilter1D(int m, double bandwidth, unsigned maxDegree, unsigned numberOfGridPoints, double xmin, double xmax, const BoundaryHandling &boundaryMethod, const unsigned char *exclusionMask=0, bool excludeCentralPoint=false)
Definition: AbsFilter1DBuilder.hh:129