npstat is hosted by Hepforge, IPPP Durham
NPStat  5.10.0
AbsFilter1DBuilder.hh
Go to the documentation of this file.
1 #ifndef NPSTAT_ABSFILTER1DBUILDER_HH_
2 #define NPSTAT_ABSFILTER1DBUILDER_HH_
3 
4 /*!
5 // \file AbsFilter1DBuilder.hh
6 //
7 // \brief Abstract interface for building local polynomial filter weights in 1-d
8 //
9 // Author: I. Volobouev
10 //
11 // November 2009
12 */
13 
14 #include <vector>
15 #include <stdexcept>
16 
17 #include "geners/ClassId.hh"
18 #include "geners/CPP11_auto_ptr.hh"
19 
21 
24 
25 namespace npstat {
26  class OrthoPoly1D;
27 
28  /**
29  // Besides providing a table of weights, the filter should remember
30  // for which point it was constructed (this info can potentially be
31  // used in subsequent cross-validation calculations)
32  */
33  class PolyFilter1D : public std::vector<long double>
34  {
35  public:
36  inline explicit PolyFilter1D(const unsigned peak) : peak_(peak) {}
37 
38  inline unsigned peakPosition() const {return peak_;}
39 
40  inline bool operator==(const PolyFilter1D& r) const
41  {return peak_ == r.peak_ &&
42  static_cast<const std::vector<long double>&>(*this) ==
43  static_cast<const std::vector<long double>&>(r);}
44 
45  inline bool operator!=(const PolyFilter1D& r) const
46  {return !(*this == r);}
47 
48  // Methods needed for I/O
49  inline gs::ClassId classId() const {return gs::ClassId(*this);}
50  bool write(std::ostream& os) const;
51 
52  static inline const char* classname() {return "npstat::PolyFilter1D";}
53  static inline unsigned version() {return 1;}
54  static PolyFilter1D* read(const gs::ClassId& id, std::istream& in);
55 
56  private:
57  PolyFilter1D();
58  unsigned peak_;
59 
60 #ifdef SWIG
61  public:
62  inline double valueAt(const unsigned i) const {return this->at(i);}
63 #endif // SWIG
64  };
65 
66  /**
67  // Abstract interface class for building local polynomial filter
68  // weights in 1-d
69  */
71  {
72  inline virtual ~AbsFilter1DBuilder() {}
73 
74  /**
75  // Length of the filter constructed at a point deeply inside
76  // the density support region
77  */
78  virtual unsigned centralWeightLength() const = 0;
79 
80  /**
81  // Should we keep all filters or can we assume that
82  // filters deeply inside the density support region
83  // are identical?
84  */
85  virtual bool keepAllFilters() const = 0;
86 
87  /**
88  // Build the filter from the given taper function, maximum
89  // polynomial degree, bin number for which this filter is
90  // constructed, and expected length of the data. The filter
91  // is constructed on the heap and later must be deleted.
92  */
93  virtual PolyFilter1D* makeFilter(const double* taper,
94  unsigned maxDegree,
95  unsigned binnum,
96  unsigned datalen) const = 0;
97 
98  /**
99  // Some filter builders may adjust boundary kernels by
100  // stretching them. In this case it may be interesting
101  // to see the bandwidth factor used.
102  */
103  inline virtual double lastBandwidthFactor() const {return 1.0;}
104 
105 #ifdef SWIG
106  inline PolyFilter1D* makeFilter_2(const std::vector<double>& taper,
107  const unsigned maxDegree,
108  const unsigned binnum,
109  const unsigned datalen) const
110  {
111  const double* tap = 0;
112  if (!taper.empty())
113  {
114  if (taper.size() < maxDegree + 1U) throw std::invalid_argument(
115  "In npstat::AbsFilter1DBuilder::makeFilter_2: "
116  "insufficient taper length");
117  tap = &taper[0];
118  }
119  return this->makeFilter(tap, maxDegree, binnum, datalen);
120  }
121 #endif // SWIG
122  };
123 
124  /**
125  // Abstract interface class for building local polynomial filter
126  // weights in 1-d via orthogonal polynomial systems
127  */
129  {
130  inline virtual ~OrthoPolyFilter1DBuilder() {}
131 
132  /** Implemented from the base */
133  virtual PolyFilter1D* makeFilter(const double* taper,
134  unsigned maxDegree,
135  unsigned binnum,
136  unsigned datalen) const;
137  /**
138  // Build the orthogonal polynomial system that can later
139  // be used to construct filters with different tapers.
140  // This is constructed on the heap and later must be deleted.
141  */
142  virtual OrthoPoly1D* makeOrthoPoly(unsigned maxDegree,
143  unsigned binnum,
144  unsigned datalen,
145  unsigned* filterCenter) const = 0;
146  };
147 
148  /**
149  // Abstract base clase for adjusting the weight function bandwidth near
150  // the boundaries by keeping the value of a certain criterion constant.
151  // Usually, this criterion is some kind of a functional (e.g., the
152  // standard deviation of the weight function) which can be changed
153  // by increasing the bandwidth.
154  */
156  {
157  public:
158  /**
159  // This class will not own the AbsDistribution1D object.
160  // It is the responsibility of the user of this class
161  // to make sure that the AbsDistribution1D object stays
162  // alive while this object is in use.
163  //
164  // It will be assumed that the distribution is symmetric
165  // about x = 0.0.
166  */
168  double centralStepSize,
169  const unsigned char* exclusionMask = 0,
170  unsigned exclusionMaskLen = 0,
171  bool excludeCentralPoint = false);
172 
173  inline virtual ~AbsBoundaryFilter1DBuilder() {}
174 
175  inline virtual unsigned centralWeightLength() const {return maxlen_;}
176 
177  inline virtual bool keepAllFilters() const
178  {return !exclusionMask_.empty();}
179 
180  virtual OrthoPoly1D* makeOrthoPoly(unsigned maxDegree,
181  unsigned binnum,
182  unsigned datalen,
183  unsigned* filterCenter) const;
184 
185  /**
186  // Return "true" for methods which fold the weight function
187  // at the boundary
188  */
189  virtual bool isFolding() const = 0;
190 
191  /**
192  // "lastBandwidthFactor" will be meaningfully defined after calling
193  // the "makeOrthoPoly" method
194  */
195  inline virtual double lastBandwidthFactor() const
196  {return lastBandwidthFactor_;}
197 
198  protected:
199  void scanTheDensity(
200  const AbsDistribution1D* distro, double h,
201  int datalen, int weightCenterPos,
202  double stepSize, double* workbuf,
203  unsigned* firstWeightUsed=0, unsigned* sizeNeeded=0) const;
204 
205  private:
207 
208  /**
209  // This function calculates the criterion which this class
210  // will attempt to keep constant throughout the whole range
211  // of x_fit by changing h
212  */
213  virtual double calculateCriterion(
214  const AbsDistribution1D* distro, double h,
215  int datalen, int weightCenterPos,
216  double stepSize, double* workbuf) const = 0;
217 
218  class Fnc;
219  friend class Fnc;
220  class Fnc : public Functor1<double,double>
221  {
222  public:
223  inline Fnc(const AbsBoundaryFilter1DBuilder& ref,
224  const unsigned datalen, const unsigned binnum)
225  : ref_(ref), datalen_(datalen), binnum_(binnum) {}
226  inline virtual ~Fnc() {}
227 
228  inline double operator()(const double& h) const
229  {
230  return ref_.calculateCriterion(ref_.distro_, h, datalen_,
231  binnum_, ref_.step_, &ref_.w_[0]);
232  }
233 
234  private:
235  const AbsBoundaryFilter1DBuilder& ref_;
236  unsigned datalen_;
237  unsigned binnum_;
238  };
239 
240  const AbsDistribution1D* distro_;
241  double step_;
242  mutable double lastBandwidthFactor_;
243  mutable double centralIntegral_;
244  unsigned maxlen_;
245  mutable unsigned centralIntegralLen_;
246  mutable std::vector<double> w_;
247  std::vector<unsigned char> exclusionMask_;
248  bool excludeCentralPoint_;
249  };
250 
251  /**
252  // A factory method for creating AbsBoundaryFilter1DBuilder objects.
253  // Parameters are as follows:
254  //
255  // boundaryHandling -- Definition of the boundary handling method.
256  //
257  // distro -- Function to use as the weight for generating
258  // LOrPE polynomials. This weight is expected
259  // to be even and to peak at 0. Note that the
260  // filter builder will not own this pointer.
261  // It is the responsibility of the user of this
262  // code to make sure that the function exists
263  // at all times when the builder is used.
264  //
265  // stepSize -- Step size (bin width) for the data grid on
266  // which density estimation will be performed.
267  //
268  // exclusionMask -- Set values of "exclusionMask" != 0 if
269  // corresponding data points have to be
270  // excluded when weights are generated.
271  // If no exclusions are necessary, just
272  // leave this array as NULL.
273  //
274  // exclusionMaskLen -- Length of the "exclusionMask" array. If
275  // it is not 0 then it must coinside with
276  // the "datalen" argument given to all future
277  // invocations of the "makeFilter" method.
278  //
279  // excludeCentralPoint -- If "true", the central point of the
280  // weight will be set to zero. This can be
281  // useful for certain types of cross validation
282  // calculations.
283  */
284  CPP11_auto_ptr<AbsBoundaryFilter1DBuilder> getBoundaryFilter1DBuilder(
285  const BoundaryHandling& boundaryHandling,
286  const AbsDistribution1D* distro, double stepSize,
287  const unsigned char* exclusionMask = 0, unsigned exclusionMaskLen = 0,
288  bool excludeCentralPoint = false);
289 
290 #ifdef SWIG
291  // Swig does not handle smart pointers properly
292  inline AbsBoundaryFilter1DBuilder* getBoundaryFilter1DBuilder_2(
293  const BoundaryHandling& boundaryMethod,
294  const AbsDistribution1D* distro, double stepSize,
295  const std::vector<int>* exclusionMask = 0,
296  bool excludeCentralPoint = false)
297  {
298  typedef CPP11_auto_ptr<AbsBoundaryFilter1DBuilder> BPtr;
299  BPtr b;
300 
301  bool have_mask = false;
302  if (exclusionMask)
303  if (!exclusionMask->empty())
304  have_mask = true;
305  if (have_mask)
306  {
307  const unsigned n = exclusionMask->size();
308  std::vector<unsigned char> mask(n, 0);
309  for (unsigned i=0; i<n; ++i)
310  if ((*exclusionMask)[i])
311  mask[i] = 1;
313  boundaryMethod, distro, stepSize,
314  &mask[0], n, excludeCentralPoint);
315  }
316  else
317  {
319  boundaryMethod, distro, stepSize,
320  (unsigned char*)0, 0U, excludeCentralPoint);
321  }
322  return b.release();
323  }
324 #endif // SWIG
325 }
326 
327 #endif // NPSTAT_ABSFILTER1DBUILDER_HH_
Interface definition for 1-d continuous statistical distributions.
API for LOrPE boundary handling methods.
Interface definitions and concrete simple functors for a variety of functor-based calculations.
Definition: AbsFilter1DBuilder.hh:156
virtual OrthoPoly1D * makeOrthoPoly(unsigned maxDegree, unsigned binnum, unsigned datalen, unsigned *filterCenter) const
virtual unsigned centralWeightLength() const
Definition: AbsFilter1DBuilder.hh:175
virtual double lastBandwidthFactor() const
Definition: AbsFilter1DBuilder.hh:195
virtual bool isFolding() const =0
AbsBoundaryFilter1DBuilder(const AbsDistribution1D *distro, double centralStepSize, const unsigned char *exclusionMask=0, unsigned exclusionMaskLen=0, bool excludeCentralPoint=false)
virtual bool keepAllFilters() const
Definition: AbsFilter1DBuilder.hh:177
Definition: BoundaryHandling.hh:21
Definition: OrthoPoly1D.hh:29
Definition: AbsFilter1DBuilder.hh:34
Definition: AbsArrayProjector.hh:14
CPP11_auto_ptr< AbsBoundaryFilter1DBuilder > getBoundaryFilter1DBuilder(const BoundaryHandling &boundaryHandling, const AbsDistribution1D *distro, double stepSize, const unsigned char *exclusionMask=0, unsigned exclusionMaskLen=0, bool excludeCentralPoint=false)
Definition: AbsDistribution1D.hh:31
Definition: AbsFilter1DBuilder.hh:71
virtual PolyFilter1D * makeFilter(const double *taper, unsigned maxDegree, unsigned binnum, unsigned datalen) const =0
virtual bool keepAllFilters() const =0
virtual unsigned centralWeightLength() const =0
virtual double lastBandwidthFactor() const
Definition: AbsFilter1DBuilder.hh:103
Definition: SimpleFunctors.hh:58
Definition: AbsFilter1DBuilder.hh:129
virtual PolyFilter1D * makeFilter(const double *taper, unsigned maxDegree, unsigned binnum, unsigned datalen) const
virtual OrthoPoly1D * makeOrthoPoly(unsigned maxDegree, unsigned binnum, unsigned datalen, unsigned *filterCenter) const =0