npstat is hosted by Hepforge, IPPP Durham
NPStat  5.10.0
ContOrthoPoly1D.hh
Go to the documentation of this file.
1 #ifndef NPSTAT_CONTORTHOPOLY1D_HH_
2 #define NPSTAT_CONTORTHOPOLY1D_HH_
3 
4 /*!
5 // \file ContOrthoPoly1D.hh
6 //
7 // \brief Continuous orthogonal polynomial series in one dimension
8 // generated for a discrete measure
9 //
10 // Author: I. Volobouev
11 //
12 // May 2017
13 */
14 
15 #include <vector>
16 #include <utility>
17 #include <cassert>
18 #include <cmath>
19 #include <map>
20 
21 #include "geners/CPP11_auto_ptr.hh"
22 
24 #include "npstat/nm/Matrix.hh"
26 #include "npstat/nm/PreciseType.hh"
28 
29 #ifdef USE_LAPACK_QUAD
31 #else
33 #endif
34 
35 namespace npstat {
36  class StorablePolySeries1D;
37  class ContOrthoPoly1DDeg;
38 
40  {
41  public:
42  typedef npstat::PreciseType<double>::type precise_type;
43 
44  // The first element of the pair is the coordinate and the
45  // second element is the weight (all weights must be non-negative)
46  typedef std::pair<double,double> MeasurePoint;
47 
48 #ifdef USE_LAPACK_QUAD
50 #else
52 #endif
54 
55  /**
56  // Main constructor, with obvious arguments. The first element
57  // of the measure pair is the coordinate and the second element
58  // is the weight (all weights must be non-negative). Internally,
59  // the measure will be sorted in the order of increasing weight
60  // and the measure coordinates will normally be shifted so that
61  // their weighted mean is at 0.
62  */
63  template<typename Numeric1, typename Numeric2>
65  const std::vector<std::pair<Numeric1,Numeric2> >& measure,
66  OrthoPolyMethod m = OPOLY_STIELTJES,
67  bool shiftMeasureCoords = true);
68 
69  /** Constructor which assumes that all weights are 1.0 */
70  template<typename Numeric>
72  const std::vector<Numeric>& coords,
73  OrthoPolyMethod m = OPOLY_STIELTJES,
74  bool shiftMeasureCoords = true);
75 
76  /** Constructor that uses the modified Chebyshev algorithm */
77  template<typename Numeric1, typename Numeric2, class OrthoPoly>
79  const std::vector<std::pair<Numeric1,Numeric2> >& measure,
80  const OrthoPoly& ops, bool shiftMeasureCoords = false);
81 
82  /**
83  // Constructor that uses the modified Chebyshev algorithm
84  // and assumes that all weights are 1.0
85  */
86 #ifndef SWIGBUG
87  // SWIG wrapping of this constructor is disabled for now.
88  // It is not obvious how to tell SWIG not to instantiate
89  // this constructor with template parameters <double, double>
90  // when the main constructor is instantiated with these parameters.
91  template<typename Numeric, class OrthoPoly>
93  const std::vector<Numeric>& coords,
94  const OrthoPoly& ops, bool shiftMeasureCoords = false);
95 #endif
96 
97  /**
98  // Function added for interface compatibility with
99  // AbsClassicalOrthoPoly1D
100  */
101  inline ContOrthoPoly1D* clone() const {return new ContOrthoPoly1D(*this);}
102 
103  //@{
104  /** A simple inspector of object properties */
105  inline unsigned maxDegree() const {return maxdeg_;}
106  inline unsigned long measureLength() const {return measure_.size();}
107  inline double weight(const unsigned long index) const
108  {return measure_.at(index).second;}
109  inline precise_type sumOfWeights() const {return wsum_;}
110  inline precise_type sumOfWeightSquares() const {return wsumsq_;}
111  inline bool areAllWeightsEqual() const {return allWeightsEqual_;}
112  inline double minCoordinate() const {return minX_;}
113  inline double maxCoordinate() const {return maxX_;}
114  inline precise_type measureShift() const {return shiftX_;}
115  inline precise_type meanCoordinate() const {return meanX_;}
116  inline precise_type location() const {return meanX_;}
117  inline precise_type scale() const {return 1.0;}
118  //@}
119 
120  //@{
121  /**
122  // Note that this returns the internal, shifted coordinate.
123  // Add the mean coordinate to restore the original coordinate.
124  */
125  inline MeasurePoint measure(const unsigned long index) const
126  {return measure_.at(index);}
127  inline double coordinate(const unsigned long index) const
128  {return measure_.at(index).first;}
129  //@}
130 
131  /** Kish's effective sample size for the measure */
132  double effectiveSampleSize() const;
133 
134  /** Return the value of one of the normalized polynomials */
135  double poly(unsigned deg, double x) const;
136 
137  /** Return the values of two of the normalized polynomials */
138  std::pair<double,double> polyPair(
139  unsigned deg1, unsigned deg2, double x) const;
140 
141  /**
142  // Return the values of all orthonormal polynomials up to some
143  // maximum degree. Size of the "polyValues" array should be
144  // at least maxdeg + 1.
145  */
146  template <typename Numeric>
147  inline void allPolys(const double x, const unsigned maxdeg,
148  Numeric* polyValues) const
149  {getAllPolys(x - shiftX_, maxdeg, polyValues);}
150 
151  /**
152  // Function added for interface compatibility with
153  // AbsClassicalOrthoPoly1D
154  */
155  inline void allpoly(const long double x, long double* values,
156  const unsigned maxdeg) const
157  {getAllPolys(x - shiftX_, maxdeg, values);}
158 
159  /**
160  // Average values of the polynomials for a sample of points.
161  // Size of the "averages" array should be at least maxdeg + 1.
162  // Note that the resulting averages are not weighted by the measure.
163  */
164  template <typename Numeric>
165  void sampleAverages(const Numeric* coords, unsigned long lenCoords,
166  double* averages, unsigned maxdeg) const;
167 
168  /**
169  // Averages of the polynomial values for the given sample of
170  // weighted points (not weighting by the measure). The first
171  // element of the pair will be treated as the coordinate
172  // and the second element will be treated as weight. Weights
173  // must not be negative. The length of array "averages" should
174  // be at least maxdeg + 1.
175  */
176  template <class Pair>
177  void weightedPointsAverages(const Pair* points, unsigned long numPoints,
178  double* averages, unsigned maxdeg) const;
179 
180  /**
181  // Average values of the pairwise polynomial products for
182  // a sample of points. The returned matrix will be symmetric
183  // and will have the dimensions (maxdeg + 1) x (maxdeg + 1).
184  // Note that the resulting averages are not weighted by the measure.
185  */
186  template <typename Numeric>
187  Matrix<double> sampleProductAverages(const Numeric* coords,
188  unsigned long lenCoords,
189  unsigned maxdeg) const;
190 
191  /** Return the value of the orthonormal polynomial series */
192  template <typename Numeric>
193  double series(const Numeric* coeffs, unsigned deg, double x) const;
194 
195  /**
196  // Build the coefficients of the orthogonal polynomial series
197  // for the given function. The length of the array "coeffs"
198  // should be at least "maxdeg" + 1. Note that the coefficients
199  // are returned in the order of increasing degree (same order
200  // is used by the "series" function).
201  */
202  template <class Functor, typename Numeric>
203  void calculateCoeffs(const Functor& fcn,
204  Numeric* coeffs, unsigned maxdeg) const;
205 
206  /**
207  // Build the coefficients of the orthogonal polynomial series
208  // for the given function in such a way that the integral of
209  // the truncated series on the [xmin, xmax] interval is constrained
210  // to "integralValue". The length of the array "coeffs" should be
211  // at least "maxdeg" + 1. Note that the coefficients are returned
212  // in the order of increasing degree (same order is used by the
213  // "series" function). The construction minimizes the ISE weighted
214  // by the empirical density.
215  //
216  // The function returns the chi-square from the corresponding
217  // constrained least squares problem. This is the sum of
218  // (coeffs[i] - c[i])^2, 0 <= i <= maxdeg, where c[i] are
219  // determined by the "calculateCoeffs" method.
220  */
221  template <class Functor>
222  double constrainedCoeffs(const Functor& fcn,
223  double *coeffs, unsigned maxdeg,
224  double xmin, double xmax,
225  double integralValue = 1.0) const;
226 
227  /**
228  // Build the coefficients of the orthogonal polynomial series
229  // for the discrete weight values (that is, fcn(x_i) = w_i,
230  // using the terminology of the "calculateCoeffs" method). This
231  // can sometimes be useful for weighted measures. Of course,
232  // for unweighted measures this results in just delta_{deg,0}.
233  */
234  void weightCoeffs(double *coeffs, unsigned maxdeg) const;
235 
236  /**
237  // Integrated squared error between the given function and the
238  // polynomial series. Parameter "integrationRulePoints" specifies
239  // the parameter "npoints" for the GaussLegendreQuadrature class.
240  // If "integrationRulePoints" is 0, the rule will be picked
241  // automatically in such a way that it integrates a polynomial
242  // of degree maxdeg*2 exactly.
243  */
244  template <class Functor>
245  double integratedSquaredError(const Functor& fcn,
246  const double *coeffs, unsigned maxdeg,
247  double xmin, double xmax,
248  unsigned integrationRulePoints=0U) const;
249 
250  /**
251  // Squared error between the given function and the polynomial
252  // series, weighted according to the measure
253  */
254  template <class Functor>
255  double weightedSquaredError(const Functor& fcn, const double *coeffs,
256  unsigned maxdeg) const;
257 
258  /**
259  // This method is useful for testing the numerical precision
260  // of the orthonormalization procedure. It returns the scalar
261  // products between various polynomials.
262  */
263  precise_type empiricalKroneckerDelta(unsigned deg1, unsigned deg2) const;
264 
265  /**
266  // A faster way to generate Kronecker deltas if the whole matrix
267  // of them is needed. The argument must not exceed the return
268  // value of the "maxDegree" function. The resulting matrix will
269  // be dimensioned (maxdeg + 1) x (maxdeg + 1).
270  */
272 
273  /**
274  // If the Kronecker deltas were calculated from a sample, they
275  // would be random and would have a covariance matrix. This
276  // is an estimate of the terms in that covariance matrix. The
277  // covariance is between delta_{deg1,deg2} and delta_{deg3,deg4}.
278  */
279  double empiricalKroneckerCovariance(unsigned deg1, unsigned deg2,
280  unsigned deg3, unsigned deg4) const;
281 
282  /**
283  // Generate principal minor of order n of the Jacobi matrix.
284  // n must not exceed "maxDegree". Note that, if you calculate
285  // roots using this Jacobi matrix, they will be shifted
286  // by meanCoordinate() (that is, add meanCoordinate() to all
287  // such roots).
288  */
290 
291  /**
292  // Integral of an orthonormal polynomials to some power without
293  // the weight (so this is not an inner product with the poly of
294  // degree 0).
295  */
296  precise_type integratePoly(unsigned degree, unsigned power,
297  double xmin, double xmax) const;
298 
299  /**
300  // Integral of the product of three orthonormal polynomials
301  // without the weight (so this is not an inner product)
302  */
303  precise_type integrateTripleProduct(unsigned deg1, unsigned deg2,
304  unsigned deg3, double xmin,
305  double xmax) const;
306 
307  /**
308  // Unweighted triple product integrals arranged in a matrix.
309  // The returned matrix will be dimensioned pairs.size() x (maxdeg+1).
310  */
312  const std::vector<std::pair<unsigned,unsigned> >& pairs,
313  unsigned maxdeg, double xmin, double xmax) const;
314 
315  /**
316  // Integral of the product of two orthonormal polynomials
317  // with some external weight function (e.g., oracle density).
318  // "EW" in the method name stands for "externally weighted".
319  //
320  // "integrationRulePoints" argument will be passed to the
321  // GaussLegendreQuadrature constructor and must be supported
322  // by that class.
323  */
324  template <class Functor>
325  double integrateEWPolyProduct(const Functor& weight,
326  unsigned deg1, unsigned deg2,
327  double xmin, double xmax,
328  unsigned integrationRulePoints) const;
329 
330  /**
331  // Integral of the polynomials with some external weight
332  // function (e.g., reconstructed or oracle density). The
333  // length of array "averages" should be at least maxdeg + 1.
334  */
335  template <class Functor>
336  void extWeightAverages(const Functor& weight,
337  const AbsIntervalQuadrature1D& quad,
338  unsigned maxdeg, double xmin, double xmax,
339  double* averages) const;
340 
341  /**
342  // Integral of all products of two orthonormal polynomials
343  // with some external weight function (e.g., oracle density).
344  */
345  template <class Functor>
347  const AbsIntervalQuadrature1D& quad,
348  unsigned maxdeg,
349  double xmin, double xmax) const;
350 
351  /**
352  // A measure-weighted average of orthonormal poly values
353  // to some power
354  */
355  double powerAverage(unsigned deg, unsigned power) const;
356 
357  /**
358  // A measure-weighted average of the product of two orthonormal
359  // poly values to some powers
360  */
361  double jointPowerAverage(unsigned deg1, unsigned power1,
362  unsigned deg2, unsigned power2) const;
363 
364  /**
365  // A measure-weighted average of the product of several
366  // orthonormal poly values
367  */
368  double jointAverage(const unsigned* degrees, unsigned nDegrees,
369  bool degreesSortedIncreasingOrder = false) const;
370 
371  //@{
372  /**
373  // A measure-weighted average that will be remembered internally
374  // so that another call to this function with the same arguments
375  // will return the answer quickly
376  */
377  double cachedJointAverage(unsigned deg1, unsigned deg2,
378  unsigned deg3, unsigned deg4) const;
379  double cachedJointAverage(unsigned deg1, unsigned deg2,
380  unsigned deg3, unsigned deg4,
381  unsigned deg5, unsigned deg6) const;
382  double cachedJointAverage(unsigned deg1, unsigned deg2,
383  unsigned deg3, unsigned deg4,
384  unsigned deg5, unsigned deg6,
385  unsigned deg7, unsigned deg8) const;
386  //@}
387 
388  /** Covariance between v_{ab} and v_{cd} */
389  double cov4(unsigned a, unsigned b, unsigned c, unsigned d) const;
390 
391  /** Covariance between v_{ab}*v_{cd} and v_{ef} */
392  double cov6(unsigned a, unsigned b, unsigned c, unsigned d,
393  unsigned e, unsigned f) const;
394 
395  /** Covariance between v_{ab}*v_{cd} and v_{ef}*v_{gh} */
396  double cov8(unsigned a, unsigned b, unsigned c, unsigned d,
397  unsigned e, unsigned f, unsigned g, unsigned h) const;
398 
399  /** Covariance between two cov4 estimates (i.e., error on the error) */
400  double covCov4(unsigned a, unsigned b, unsigned c, unsigned d,
401  unsigned e, unsigned f, unsigned g, unsigned h) const;
402 
403  /**
404  // Alternative implementation of "cov8". It is slower than "cov8"
405  // but easier to program and verify. Useful mainly for testing "cov8".
406  */
407  double slowCov8(unsigned a, unsigned b, unsigned c, unsigned d,
408  unsigned e, unsigned f, unsigned g, unsigned h) const;
409 
410  /** Estimate of eps_{mn} expectation */
411  double epsExpectation(unsigned m, unsigned n, bool highOrder) const;
412 
413  /**
414  // Estimate of covariance between eps_{ij} and eps_{mn}.
415  //
416  // The result should be identical with "empiricalKroneckerCovariance"
417  // in case "highOrder" argument is "false". However, this method
418  // caches results of various calculations of averages and should
419  // perform better inside cycles over indices.
420  */
421  double epsCovariance(unsigned i, unsigned j,
422  unsigned m, unsigned n, bool highOrder) const;
423 
424  /**
425  // Covariance matrix created by multiple calls to "epsCovariance"
426  */
428  const std::vector<std::pair<unsigned,unsigned> >& pairs,
429  bool highOrder) const;
430 
431  /** Construct a StorablePolySeries1D object out of this one */
432  CPP11_auto_ptr<StorablePolySeries1D> makeStorablePolySeries(
433  double xmin, double xmax,
434  const double *coeffs=0, unsigned maxdeg=0) const;
435 
436  /** Examine recurrence coefficients */
437  inline std::pair<precise_type,precise_type>
438  recurrenceCoeffs(const unsigned k) const
439  {
440  const Recurrence& rc(rCoeffs_.at(k));
441  return std::pair<precise_type,precise_type>(rc.alpha, rc.sqrbeta);
442  }
443 
444  /**
445  // Examine recurrence coefficients for the corresponding
446  // monic polynomials
447  */
448  inline std::pair<precise_type,precise_type>
449  monicRecurrenceCoeffs(const unsigned k) const
450  {
451  const Recurrence& rc(rCoeffs_.at(k));
452  return std::pair<precise_type,precise_type>(rc.alpha, rc.beta);
453  }
454 
455  private:
456  static const precise_type precise_zero;
457 
458  struct Recurrence
459  {
460  inline Recurrence(const precise_type a,
461  const precise_type b)
462  : alpha(a), beta(b)
463  {
464  assert(beta > 0.0);
465  sqrbeta = std::sqrt(beta);
466  }
467 
468  precise_type alpha;
469  precise_type beta;
470  precise_type sqrbeta;
471  };
472 
473  class PolyFcn;
474  friend class PolyFcn;
475 
476  class PolyFcn : public Functor1<precise_type, precise_type>
477  {
478  public:
479  inline PolyFcn(const ContOrthoPoly1D& p,
480  const unsigned d1, const unsigned power)
481  : poly(p), deg1(d1), polypow(power) {}
482 
483  inline virtual ~PolyFcn() {}
484 
485  inline precise_type operator()(const precise_type& x) const
486  {
487  precise_type p = 1.0;
488  if (polypow)
489  {
490  p = poly.normpoly(deg1, x);
491  switch (polypow)
492  {
493  case 1U:
494  break;
495  case 2U:
496  p *= p;
497  break;
498  default:
499  const precise_type myPow = polypow;
500  p = std::pow(p, myPow);
501  break;
502  }
503  }
504  return p;
505  }
506 
507  private:
508  const ContOrthoPoly1D& poly;
509  unsigned deg1;
510  unsigned polypow;
511  };
512 
513  class TripleProdFcn;
514  friend class TripleProdFcn;
515 
516  class TripleProdFcn : public Functor1<precise_type, precise_type>
517  {
518  public:
519  inline TripleProdFcn(const ContOrthoPoly1D& p, const unsigned d1,
520  const unsigned d2, const unsigned d3)
521  : poly(p)
522  {
523  degs[0] = d1;
524  degs[1] = d2;
525  degs[2] = d3;
526  degmax = std::max(d1, std::max(d2, d3));
527  }
528 
529  inline virtual ~TripleProdFcn() {}
530 
531  inline precise_type operator()(const precise_type& x) const
532  {return poly.normpolyprod(degs, 3, degmax, x);}
533 
534  private:
535  const ContOrthoPoly1D& poly;
536  unsigned degs[3];
537  unsigned degmax;
538  };
539 
540  template<class Funct> class EWPolyProductFcn;
541  template<class Funct> friend class EWPolyProductFcn;
542 
543  template<class Funct>
544  class EWPolyProductFcn : public Functor1<precise_type, precise_type>
545  {
546  public:
547  inline EWPolyProductFcn(const ContOrthoPoly1D& p, const Funct& w,
548  const unsigned d1, const unsigned d2)
549  : poly(p), wf(w)
550  {
551  degs[0] = d1;
552  degs[1] = d2;
553  degmax = std::max(d1, d2);
554  }
555 
556  inline virtual ~EWPolyProductFcn() {}
557 
558  inline precise_type operator()(const precise_type& x) const
559  {return poly.normpolyprod(degs,2,degmax,x-poly.shiftX_)*wf(x);}
560 
561  private:
562  const ContOrthoPoly1D& poly;
563  const Funct& wf;
564  unsigned degs[2];
565  unsigned degmax;
566  };
567 
568  class MemoKey
569  {
570  public:
571  inline MemoKey(const unsigned d0, const unsigned d1,
572  const unsigned d2, const unsigned d3)
573  : nDegs_(4U)
574  {
575  degs_[0] = d0;
576  degs_[1] = d1;
577  degs_[2] = d2;
578  degs_[3] = d3;
579  std::sort(degs_, degs_+nDegs_);
580  }
581 
582  inline MemoKey(const unsigned d0, const unsigned d1,
583  const unsigned d2, const unsigned d3,
584  const unsigned d4, const unsigned d5)
585  : nDegs_(6U)
586  {
587  degs_[0] = d0;
588  degs_[1] = d1;
589  degs_[2] = d2;
590  degs_[3] = d3;
591  degs_[4] = d4;
592  degs_[5] = d5;
593  std::sort(degs_, degs_+nDegs_);
594  }
595 
596  inline MemoKey(const unsigned d0, const unsigned d1,
597  const unsigned d2, const unsigned d3,
598  const unsigned d4, const unsigned d5,
599  const unsigned d6, const unsigned d7)
600  : nDegs_(8U)
601  {
602  degs_[0] = d0;
603  degs_[1] = d1;
604  degs_[2] = d2;
605  degs_[3] = d3;
606  degs_[4] = d4;
607  degs_[5] = d5;
608  degs_[6] = d6;
609  degs_[7] = d7;
610  std::sort(degs_, degs_+nDegs_);
611  }
612 
613  inline const unsigned* degrees() const {return degs_;}
614  inline unsigned nDegrees() const {return nDegs_;}
615 
616  inline bool operator<(const MemoKey& r) const
617  {
618  if (nDegs_ < r.nDegs_)
619  return true;
620  if (r.nDegs_ < nDegs_)
621  return false;
622  for (unsigned i=0; i<nDegs_; ++i)
623  {
624  if (degs_[i] < r.degs_[i])
625  return true;
626  if (r.degs_[i] < degs_[i])
627  return false;
628  }
629  return false;
630  }
631 
632  inline bool operator==(const MemoKey& r) const
633  {
634  if (nDegs_ != r.nDegs_)
635  return false;
636  for (unsigned i=0; i<nDegs_; ++i)
637  if (degs_[i] != r.degs_[i])
638  return false;
639  return true;
640  }
641 
642  inline bool operator!=(const MemoKey& r) const
643  {return !(*this == r);}
644 
645  private:
646  unsigned degs_[8];
647  unsigned nDegs_;
648  };
649 
650  template<typename Numeric1, typename Numeric2>
651  void copyMeasure(const std::vector<std::pair<Numeric1,Numeric2> >& inMeasure);
652 
653  template<typename Numeric>
654  void copyInputCoords(const std::vector<Numeric>& inCoords);
655 
656  void calcMeanXandWeightSums(bool shiftTheMeasure);
657  void calcRecurrenceCoeffs(OrthoPolyMethod m);
658 
659  void calcRecurrenceStieltjes();
660  void calcRecurrenceLanczos();
661 
662  template<class OrthoPoly>
663  long double sigma(const OrthoPoly& ops, int k, int l,
664  Matrix<long double>& sigMatrix,
665  Matrix<int>& flagMatrix) const;
666 
667  template<class OrthoPoly>
668  void calculateMonicAverages(const OrthoPoly& ops, long double* averages,
669  unsigned maxdeg);
670 
671  template<class OrthoPoly>
672  void calcRecurrenceModifiedChebyshev(const OrthoPoly& ops);
673 
674  precise_type monicpoly(unsigned degree, precise_type x) const;
675 
676  std::pair<precise_type,precise_type> monicInnerProducts(
677  unsigned degree) const;
678 
679  precise_type normpoly(unsigned degree, precise_type x) const;
680 
681  std::pair<precise_type,precise_type> twonormpoly(
682  unsigned deg1, unsigned deg2, precise_type x) const;
683 
684  precise_type normpolyprod(const unsigned* degrees, unsigned nDeg,
685  unsigned degmax, precise_type x) const;
686 
687  template <typename Numeric>
688  precise_type normseries(const Numeric* coeffs, unsigned maxdeg,
689  precise_type x) const;
690 
691  template <typename Numeric>
692  void getAllPolys(double x, unsigned maxdeg, Numeric *polyValues) const;
693 
694  double getCachedAverage(const MemoKey& key) const;
695 
696  std::vector<MeasurePoint> measure_;
697  std::vector<Recurrence> rCoeffs_;
698  mutable std::map<MemoKey,double> cachedAverages_;
699  precise_type wsum_;
700  precise_type wsumsq_;
701  precise_type meanX_;
702  precise_type shiftX_;
703  double minX_;
704  double maxX_;
705  unsigned maxdeg_;
706  bool allWeightsEqual_;
707 
708 #ifdef SWIG
709  public:
710  inline StorablePolySeries1D* makeStorablePolySeries_2(
711  double xmin, double xmax, const std::vector<double>& coeffs) const
712  {
713  CPP11_auto_ptr<StorablePolySeries1D> ptr;
714  if (coeffs.empty())
715  ptr = this->makeStorablePolySeries(xmin, xmax);
716  else
717  ptr = this->makeStorablePolySeries(xmin, xmax, &coeffs[0],
718  coeffs.size() - 1);
719  return ptr.release();
720  }
721 
722  inline double sumOfWeights_2() const
723  {return sumOfWeights();}
724 
725  inline double sumOfWeightSquares_2() const
726  {return sumOfWeightSquares();}
727 
728  inline double meanCoordinate_2() const
729  {return meanCoordinate();}
730 
731  inline double empiricalKroneckerDelta_2(const unsigned deg1,
732  const unsigned deg2) const
733  {return empiricalKroneckerDelta(deg1, deg2);}
734 
735  inline Matrix<double> jacobiMatrix_2(const unsigned n) const
736  {return jacobiMatrix(n);}
737 
738  inline double integratePoly_2(const unsigned degree, const unsigned power,
739  const double xmin, const double xmax) const
740  {return integratePoly(degree, power, xmin, xmax);}
741 
742  inline double integrateTripleProduct_2(const unsigned deg1, const unsigned deg2,
743  const unsigned deg3, const double xmin,
744  const double xmax) const
745  {return integrateTripleProduct(deg1, deg2, deg3, xmin, xmax);}
746 
747  inline Matrix<double> tripleProductMatrix_2(
748  const std::vector<std::pair<unsigned,unsigned> >& pairs,
749  const unsigned maxdeg, const double xmin, const double xmax) const
750  {return tripleProductMatrix(pairs, maxdeg, xmin, xmax);}
751 
752  inline double location2() const {return location();}
753 
754  inline double scale2() const {return scale();}
755 
756  inline std::pair<double,double> monicRecurrenceCoeffs2(const unsigned deg) const
757  {
758  const std::pair<precise_type,precise_type>& p = monicRecurrenceCoeffs(deg);
759  return std::pair<double,double>(p.first, p.second);
760  }
761 
762  inline std::pair<double,double> recurrenceCoeffs2(const unsigned deg) const
763  {
764  const std::pair<precise_type,precise_type>& p = recurrenceCoeffs(deg);
765  return std::pair<double,double>(p.first, p.second);
766  }
767 #endif // SWIG
768  };
769 
770  /**
771  // A functor for a particular degree polynomial of the given
772  // ContOrthoPoly1D. The poly system is not copied, only a reference
773  // is used. It is a responsibility of the user to make sure that the
774  // lifetime of the poly system object exceeds the lifetime of the functor.
775  */
776  class ContOrthoPoly1DDeg : public Functor1<long double, long double>
777  {
778  public:
779  inline ContOrthoPoly1DDeg(const ContOrthoPoly1D& fcn,
780  const unsigned degree,
781  const long double normfactor=1.0L)
782  : fcn_(fcn), norm_(normfactor), deg_(degree)
783  {
784  if (deg_ > fcn_.maxDegree()) throw std::invalid_argument(
785  "In npstat::ContOrthoPoly1DDeg::constructor: "
786  "degree argument is out of range");
787  }
788 
789  inline virtual ~ContOrthoPoly1DDeg() {}
790 
791  inline virtual long double operator()(const long double& a) const
792  {return norm_*fcn_.poly(deg_, a);}
793 
794  private:
796 
797  const ContOrthoPoly1D& fcn_;
798  long double norm_;
799  unsigned deg_;
800  };
801 }
802 
803 #include "npstat/nm/ContOrthoPoly1D.icc"
804 
805 #endif // NPSTAT_CONTORTHOPOLY1D_HH_
Base class for quadratures on the [-1, 1] interval.
bool operator==(const npstat::BoxND< Numeric > &l, const npstat::BoxND< Numeric > &r)
Gauss-Legendre quadratures in quad precision.
Gauss-Legendre quadratures in long double precision.
Template matrix class.
Enumeration of methods used to create orthogonal polynomials with discrete weights.
Compile-time deduction of an appropriate precise numeric type.
Interface definitions and concrete simple functors for a variety of functor-based calculations.
Definition: AbsIntervalQuadrature1D.hh:19
Definition: ContOrthoPoly1D.hh:777
Definition: ContOrthoPoly1D.hh:40
CPP11_auto_ptr< StorablePolySeries1D > makeStorablePolySeries(double xmin, double xmax, const double *coeffs=0, unsigned maxdeg=0) const
MeasurePoint measure(const unsigned long index) const
Definition: ContOrthoPoly1D.hh:125
ContOrthoPoly1D(unsigned maxDegree, const std::vector< Numeric > &coords, OrthoPolyMethod m=OPOLY_STIELTJES, bool shiftMeasureCoords=true)
ContOrthoPoly1D(unsigned maxDegree, const std::vector< Numeric > &coords, const OrthoPoly &ops, bool shiftMeasureCoords=false)
double integrateEWPolyProduct(const Functor &weight, unsigned deg1, unsigned deg2, double xmin, double xmax, unsigned integrationRulePoints) const
void sampleAverages(const Numeric *coords, unsigned long lenCoords, double *averages, unsigned maxdeg) const
double empiricalKroneckerCovariance(unsigned deg1, unsigned deg2, unsigned deg3, unsigned deg4) const
double cov8(unsigned a, unsigned b, unsigned c, unsigned d, unsigned e, unsigned f, unsigned g, unsigned h) const
Matrix< double > epsCovarianceMatrix(const std::vector< std::pair< unsigned, unsigned > > &pairs, bool highOrder) const
double epsExpectation(unsigned m, unsigned n, bool highOrder) const
double integratedSquaredError(const Functor &fcn, const double *coeffs, unsigned maxdeg, double xmin, double xmax, unsigned integrationRulePoints=0U) const
void weightedPointsAverages(const Pair *points, unsigned long numPoints, double *averages, unsigned maxdeg) const
double slowCov8(unsigned a, unsigned b, unsigned c, unsigned d, unsigned e, unsigned f, unsigned g, unsigned h) const
ContOrthoPoly1D(unsigned maxDegree, const std::vector< std::pair< Numeric1, Numeric2 > > &measure, const OrthoPoly &ops, bool shiftMeasureCoords=false)
Matrix< double > sampleProductAverages(const Numeric *coords, unsigned long lenCoords, unsigned maxdeg) const
void extWeightAverages(const Functor &weight, const AbsIntervalQuadrature1D &quad, unsigned maxdeg, double xmin, double xmax, double *averages) const
double covCov4(unsigned a, unsigned b, unsigned c, unsigned d, unsigned e, unsigned f, unsigned g, unsigned h) const
double constrainedCoeffs(const Functor &fcn, double *coeffs, unsigned maxdeg, double xmin, double xmax, double integralValue=1.0) const
ContOrthoPoly1D(unsigned maxDegree, const std::vector< std::pair< Numeric1, Numeric2 > > &measure, OrthoPolyMethod m=OPOLY_STIELTJES, bool shiftMeasureCoords=true)
void allPolys(const double x, const unsigned maxdeg, Numeric *polyValues) const
Definition: ContOrthoPoly1D.hh:147
precise_type integratePoly(unsigned degree, unsigned power, double xmin, double xmax) const
double powerAverage(unsigned deg, unsigned power) const
double series(const Numeric *coeffs, unsigned deg, double x) const
std::pair< precise_type, precise_type > monicRecurrenceCoeffs(const unsigned k) const
Definition: ContOrthoPoly1D.hh:449
void weightCoeffs(double *coeffs, unsigned maxdeg) const
void allpoly(const long double x, long double *values, const unsigned maxdeg) const
Definition: ContOrthoPoly1D.hh:155
double cov4(unsigned a, unsigned b, unsigned c, unsigned d) const
Matrix< double > extWeightProductAverages(const Functor &weight, const AbsIntervalQuadrature1D &quad, unsigned maxdeg, double xmin, double xmax) const
std::pair< double, double > polyPair(unsigned deg1, unsigned deg2, double x) const
precise_type empiricalKroneckerDelta(unsigned deg1, unsigned deg2) const
double cachedJointAverage(unsigned deg1, unsigned deg2, unsigned deg3, unsigned deg4) const
Matrix< precise_type > tripleProductMatrix(const std::vector< std::pair< unsigned, unsigned > > &pairs, unsigned maxdeg, double xmin, double xmax) const
double poly(unsigned deg, double x) const
ContOrthoPoly1D * clone() const
Definition: ContOrthoPoly1D.hh:101
std::pair< precise_type, precise_type > recurrenceCoeffs(const unsigned k) const
Definition: ContOrthoPoly1D.hh:438
double weightedSquaredError(const Functor &fcn, const double *coeffs, unsigned maxdeg) const
double epsCovariance(unsigned i, unsigned j, unsigned m, unsigned n, bool highOrder) const
double cov6(unsigned a, unsigned b, unsigned c, unsigned d, unsigned e, unsigned f) const
unsigned maxDegree() const
Definition: ContOrthoPoly1D.hh:105
precise_type integrateTripleProduct(unsigned deg1, unsigned deg2, unsigned deg3, double xmin, double xmax) const
void calculateCoeffs(const Functor &fcn, Numeric *coeffs, unsigned maxdeg) const
Matrix< precise_type > jacobiMatrix(unsigned n) const
double jointAverage(const unsigned *degrees, unsigned nDegrees, bool degreesSortedIncreasingOrder=false) const
double effectiveSampleSize() const
double jointPowerAverage(unsigned deg1, unsigned power1, unsigned deg2, unsigned power2) const
Matrix< precise_type > empiricalKroneckerMatrix(unsigned maxdeg) const
Definition: GaussLegendreQuadratureQ.hh:26
Definition: GaussLegendreQuadrature.hh:27
Definition: AbsArrayProjector.hh:14
OrthoPolyMethod
Definition: OrthoPolyMethod.hh:20
Definition: SimpleFunctors.hh:58