ContOrthoPoly1D.hh
No OneTemporary
Actions

Size

31 KB

Referenced Files

None

Subscribers

None

ContOrthoPoly1D.hh
View Options

	#ifndef NPSTAT_CONTORTHOPOLY1D_HH_
	#define NPSTAT_CONTORTHOPOLY1D_HH_

	/*!
	// \file ContOrthoPoly1D.hh
	//
	// \brief Continuous orthogonal polynomial series in one dimension
	// generated for a discrete measure
	//
	// Author: I. Volobouev
	//
	// May 2017
	*/

	#include <vector>
	#include <utility>
	#include <cassert>
	#include <cmath>
	#include <map>

	#include "geners/CPP11_auto_ptr.hh"

	#include "npstat/nm/OrthoPolyMethod.hh"
	#include "npstat/nm/Matrix.hh"
	#include "npstat/nm/SimpleFunctors.hh"
	#include "npstat/nm/PreciseType.hh"
	#include "npstat/nm/AbsIntervalQuadrature1D.hh"

	#ifdef USE_LAPACK_QUAD
	#include "npstat/nm/GaussLegendreQuadratureQ.hh"
	#else
	#include "npstat/nm/GaussLegendreQuadrature.hh"
	#endif

	namespace npstat {
	class StorablePolySeries1D;
	class ContOrthoPoly1DDeg;

	class ContOrthoPoly1D
	{
	public:
	typedef npstat::PreciseType<double>::type precise_type;

	// The first element of the pair is the coordinate and the
	// second element is the weight (all weights must be non-negative)
	typedef std::pair<double,double> MeasurePoint;

	#ifdef USE_LAPACK_QUAD
	typedef GaussLegendreQuadratureQ PreciseQuadrature;
	#else
	typedef GaussLegendreQuadrature PreciseQuadrature;
	#endif
	typedef ContOrthoPoly1DDeg degree_functor;

	/**
	// Main constructor, with obvious arguments. The first element
	// of the measure pair is the coordinate and the second element
	// is the weight (all weights must be non-negative). Internally,
	// the measure will be sorted in the order of increasing weight
	// and the measure coordinates will normally be shifted so that
	// their weighted mean is at 0.
	*/
	template<typename Numeric1, typename Numeric2>
	ContOrthoPoly1D(unsigned maxDegree,
	const std::vector<std::pair<Numeric1,Numeric2> >& measure,
	OrthoPolyMethod m = OPOLY_STIELTJES,
	bool shiftMeasureCoords = true);

	/** Constructor which assumes that all weights are 1.0 */
	template<typename Numeric>
	ContOrthoPoly1D(unsigned maxDegree,
	const std::vector<Numeric>& coords,
	OrthoPolyMethod m = OPOLY_STIELTJES,
	bool shiftMeasureCoords = true);

	/** Constructor that uses the modified Chebyshev algorithm */
	template<typename Numeric1, typename Numeric2, class OrthoPoly>
	ContOrthoPoly1D(unsigned maxDegree,
	const std::vector<std::pair<Numeric1,Numeric2> >& measure,
	const OrthoPoly& ops, bool shiftMeasureCoords = false);

	/**
	// Constructor that uses the modified Chebyshev algorithm
	// and assumes that all weights are 1.0
	*/
	#ifndef SWIGBUG
	// SWIG wrapping of this constructor is disabled for now.
	// It is not obvious how to tell SWIG not to instantiate
	// this constructor with template parameters <double, double>
	// when the main constructor is instantiated with these parameters.
	template<typename Numeric, class OrthoPoly>
	ContOrthoPoly1D(unsigned maxDegree,
	const std::vector<Numeric>& coords,
	const OrthoPoly& ops, bool shiftMeasureCoords = false);
	#endif

	/**
	// Function added for interface compatibility with
	// AbsClassicalOrthoPoly1D
	*/
	inline ContOrthoPoly1D* clone() const {return new ContOrthoPoly1D(*this);}

	//@{
	/** A simple inspector of object properties */
	inline unsigned maxDegree() const {return maxdeg_;}
	inline unsigned long measureLength() const {return measure_.size();}
	inline double weight(const unsigned long index) const
	{return measure_.at(index).second;}
	inline precise_type sumOfWeights() const {return wsum_;}
	inline precise_type sumOfWeightSquares() const {return wsumsq_;}
	inline bool areAllWeightsEqual() const {return allWeightsEqual_;}
	inline double minCoordinate() const {return minX_;}
	inline double maxCoordinate() const {return maxX_;}
	inline precise_type measureShift() const {return shiftX_;}
	inline precise_type meanCoordinate() const {return meanX_;}
	inline precise_type location() const {return meanX_;}
	inline precise_type scale() const {return 1.0;}
	//@}

	//@{
	/**
	// Note that this returns the internal, shifted coordinate.
	// Add the mean coordinate to restore the original coordinate.
	*/
	inline MeasurePoint measure(const unsigned long index) const
	{return measure_.at(index);}
	inline double coordinate(const unsigned long index) const
	{return measure_.at(index).first;}
	//@}

	/** Kish's effective sample size for the measure */
	double effectiveSampleSize() const;

	/** Return the value of one of the normalized polynomials */
	double poly(unsigned deg, double x) const;

	/** Return the values of two of the normalized polynomials */
	std::pair<double,double> polyPair(
	unsigned deg1, unsigned deg2, double x) const;

	/**
	// Return the values of all orthonormal polynomials up to some
	// maximum degree. Size of the "polyValues" array should be
	// at least maxdeg + 1.
	*/
	template <typename Numeric>
	inline void allPolys(const double x, const unsigned maxdeg,
	Numeric* polyValues) const
	{getAllPolys(x - shiftX_, maxdeg, polyValues);}

	/**
	// Function added for interface compatibility with
	// AbsClassicalOrthoPoly1D
	*/
	inline void allpoly(const long double x, long double* values,
	const unsigned maxdeg) const
	{getAllPolys(x - shiftX_, maxdeg, values);}

	/**
	// Average values of the polynomials for a sample of points.
	// Size of the "averages" array should be at least maxdeg + 1.
	// Note that the resulting averages are not weighted by the measure.
	*/
	template <typename Numeric>
	void sampleAverages(const Numeric* coords, unsigned long lenCoords,
	double* averages, unsigned maxdeg) const;

	/**
	// Averages of the polynomial values for the given sample of
	// weighted points (not weighting by the measure). The first
	// element of the pair will be treated as the coordinate
	// and the second element will be treated as weight. Weights
	// must not be negative. The length of array "averages" should
	// be at least maxdeg + 1.
	*/
	template <class Pair>
	void weightedPointsAverages(const Pair* points, unsigned long numPoints,
	double* averages, unsigned maxdeg) const;

	/**
	// Average values of the pairwise polynomial products for
	// a sample of points. The returned matrix will be symmetric
	// and will have the dimensions (maxdeg + 1) x (maxdeg + 1).
	// Note that the resulting averages are not weighted by the measure.
	*/
	template <typename Numeric>
	Matrix<double> sampleProductAverages(const Numeric* coords,
	unsigned long lenCoords,
	unsigned maxdeg) const;

	/** Return the value of the orthonormal polynomial series */
	template <typename Numeric>
	double series(const Numeric* coeffs, unsigned deg, double x) const;

	/**
	// Build the coefficients of the orthogonal polynomial series
	// for the given function. The length of the array "coeffs"
	// should be at least "maxdeg" + 1. Note that the coefficients
	// are returned in the order of increasing degree (same order
	// is used by the "series" function).
	*/
	template <class Functor, typename Numeric>
	void calculateCoeffs(const Functor& fcn,
	Numeric* coeffs, unsigned maxdeg) const;

	/**
	// Build the coefficients of the orthogonal polynomial series
	// for the given function in such a way that the integral of
	// the truncated series on the [xmin, xmax] interval is constrained
	// to "integralValue". The length of the array "coeffs" should be
	// at least "maxdeg" + 1. Note that the coefficients are returned
	// in the order of increasing degree (same order is used by the
	// "series" function). The construction minimizes the ISE weighted
	// by the empirical density.
	//
	// The function returns the chi-square from the corresponding
	// constrained least squares problem. This is the sum of
	// (coeffs[i] - c[i])^2, 0 <= i <= maxdeg, where c[i] are
	// determined by the "calculateCoeffs" method.
	*/
	template <class Functor>
	double constrainedCoeffs(const Functor& fcn,
	double *coeffs, unsigned maxdeg,
	double xmin, double xmax,
	double integralValue = 1.0) const;

	/**
	// Build the coefficients of the orthogonal polynomial series
	// for the discrete weight values (that is, fcn(x_i) = w_i,
	// using the terminology of the "calculateCoeffs" method). This
	// can sometimes be useful for weighted measures. Of course,
	// for unweighted measures this results in just delta_{deg,0}.
	*/
	void weightCoeffs(double *coeffs, unsigned maxdeg) const;

	/**
	// Integrated squared error between the given function and the
	// polynomial series. Parameter "integrationRulePoints" specifies
	// the parameter "npoints" for the GaussLegendreQuadrature class.
	// If "integrationRulePoints" is 0, the rule will be picked
	// automatically in such a way that it integrates a polynomial
	// of degree maxdeg*2 exactly.
	*/
	template <class Functor>
	double integratedSquaredError(const Functor& fcn,
	const double *coeffs, unsigned maxdeg,
	double xmin, double xmax,
	unsigned integrationRulePoints=0U) const;

	/**
	// Squared error between the given function and the polynomial
	// series, weighted according to the measure
	*/
	template <class Functor>
	double weightedSquaredError(const Functor& fcn, const double *coeffs,
	unsigned maxdeg) const;

	/**
	// This method is useful for testing the numerical precision
	// of the orthonormalization procedure. It returns the scalar
	// products between various polynomials.
	*/
	precise_type empiricalKroneckerDelta(unsigned deg1, unsigned deg2) const;

	/**
	// A faster way to generate Kronecker deltas if the whole matrix
	// of them is needed. The argument must not exceed the return
	// value of the "maxDegree" function. The resulting matrix will
	// be dimensioned (maxdeg + 1) x (maxdeg + 1).
	*/
	Matrix<precise_type> empiricalKroneckerMatrix(unsigned maxdeg) const;

	/**
	// If the Kronecker deltas were calculated from a sample, they
	// would be random and would have a covariance matrix. This
	// is an estimate of the terms in that covariance matrix. The
	// covariance is between delta_{deg1,deg2} and delta_{deg3,deg4}.
	*/
	double empiricalKroneckerCovariance(unsigned deg1, unsigned deg2,
	unsigned deg3, unsigned deg4) const;

	/**
	// Generate principal minor of order n of the Jacobi matrix.
	// n must not exceed "maxDegree". Note that, if you calculate
	// roots using this Jacobi matrix, they will be shifted
	// by meanCoordinate() (that is, add meanCoordinate() to all
	// such roots).
	*/
	Matrix<precise_type> jacobiMatrix(unsigned n) const;

	/**
	// Integral of an orthonormal polynomials to some power without
	// the weight (so this is not an inner product with the poly of
	// degree 0).
	*/
	precise_type integratePoly(unsigned degree, unsigned power,
	double xmin, double xmax) const;

	/**
	// Integral of the product of three orthonormal polynomials
	// without the weight (so this is not an inner product)
	*/
	precise_type integrateTripleProduct(unsigned deg1, unsigned deg2,
	unsigned deg3, double xmin,
	double xmax) const;

	/**
	// Unweighted triple product integrals arranged in a matrix.
	// The returned matrix will be dimensioned pairs.size() x (maxdeg+1).
	*/
	Matrix<precise_type> tripleProductMatrix(
	const std::vector<std::pair<unsigned,unsigned> >& pairs,
	unsigned maxdeg, double xmin, double xmax) const;

	/**
	// Integral of the product of two orthonormal polynomials
	// with some external weight function (e.g., oracle density).
	// "EW" in the method name stands for "externally weighted".
	//
	// "integrationRulePoints" argument will be passed to the
	// GaussLegendreQuadrature constructor and must be supported
	// by that class.
	*/
	template <class Functor>
	double integrateEWPolyProduct(const Functor& weight,
	unsigned deg1, unsigned deg2,
	double xmin, double xmax,
	unsigned integrationRulePoints) const;

	/**
	// Integral of the polynomials with some external weight
	// function (e.g., reconstructed or oracle density). The
	// length of array "averages" should be at least maxdeg + 1.
	*/
	template <class Functor>
	void extWeightAverages(const Functor& weight,
	const AbsIntervalQuadrature1D& quad,
	unsigned maxdeg, double xmin, double xmax,
	double* averages) const;

	/**
	// Integral of all products of two orthonormal polynomials
	// with some external weight function (e.g., oracle density).
	*/
	template <class Functor>
	Matrix<double> extWeightProductAverages(const Functor& weight,
	const AbsIntervalQuadrature1D& quad,
	unsigned maxdeg,
	double xmin, double xmax) const;

	/**
	// A measure-weighted average of orthonormal poly values
	// to some power
	*/
	double powerAverage(unsigned deg, unsigned power) const;

	/**
	// A measure-weighted average of the product of two orthonormal
	// poly values to some powers
	*/
	double jointPowerAverage(unsigned deg1, unsigned power1,
	unsigned deg2, unsigned power2) const;

	/**
	// A measure-weighted average of the product of several
	// orthonormal poly values
	*/
	double jointAverage(const unsigned* degrees, unsigned nDegrees,
	bool degreesSortedIncreasingOrder = false) const;

	//@{
	/**
	// A measure-weighted average that will be remembered internally
	// so that another call to this function with the same arguments
	// will return the answer quickly
	*/
	double cachedJointAverage(unsigned deg1, unsigned deg2,
	unsigned deg3, unsigned deg4) const;
	double cachedJointAverage(unsigned deg1, unsigned deg2,
	unsigned deg3, unsigned deg4,
	unsigned deg5, unsigned deg6) const;
	double cachedJointAverage(unsigned deg1, unsigned deg2,
	unsigned deg3, unsigned deg4,
	unsigned deg5, unsigned deg6,
	unsigned deg7, unsigned deg8) const;
	//@}

	/** Covariance between v_{ab} and v_{cd} */
	double cov4(unsigned a, unsigned b, unsigned c, unsigned d) const;

	/** Covariance between v_{ab}v_{cd} and v_{ef} /
	double cov6(unsigned a, unsigned b, unsigned c, unsigned d,
	unsigned e, unsigned f) const;

	/** Covariance between v_{ab}v_{cd} and v_{ef}v_{gh} */
	double cov8(unsigned a, unsigned b, unsigned c, unsigned d,
	unsigned e, unsigned f, unsigned g, unsigned h) const;

	/** Covariance between two cov4 estimates (i.e., error on the error) */
	double covCov4(unsigned a, unsigned b, unsigned c, unsigned d,
	unsigned e, unsigned f, unsigned g, unsigned h) const;

	/**
	// Alternative implementation of "cov8". It is slower than "cov8"
	// but easier to program and verify. Useful mainly for testing "cov8".
	*/
	double slowCov8(unsigned a, unsigned b, unsigned c, unsigned d,
	unsigned e, unsigned f, unsigned g, unsigned h) const;

	/** Estimate of eps_{mn} expectation */
	double epsExpectation(unsigned m, unsigned n, bool highOrder) const;

	/**
	// Estimate of covariance between eps_{ij} and eps_{mn}.
	//
	// The result should be identical with "empiricalKroneckerCovariance"
	// in case "highOrder" argument is "false". However, this method
	// caches results of various calculations of averages and should
	// perform better inside cycles over indices.
	*/
	double epsCovariance(unsigned i, unsigned j,
	unsigned m, unsigned n, bool highOrder) const;

	/**
	// Covariance matrix created by multiple calls to "epsCovariance"
	*/
	Matrix<double> epsCovarianceMatrix(
	const std::vector<std::pair<unsigned,unsigned> >& pairs,
	bool highOrder) const;

	/** Construct a StorablePolySeries1D object out of this one */
	CPP11_auto_ptr<StorablePolySeries1D> makeStorablePolySeries(
	double xmin, double xmax,
	const double *coeffs=0, unsigned maxdeg=0) const;

	/** Examine recurrence coefficients */
	inline std::pair<precise_type,precise_type>
	recurrenceCoeffs(const unsigned k) const
	{
	const Recurrence& rc(rCoeffs_.at(k));
	return std::pair<precise_type,precise_type>(rc.alpha, rc.sqrbeta);
	}

	/**
	// Examine recurrence coefficients for the corresponding
	// monic polynomials
	*/
	inline std::pair<precise_type,precise_type>
	monicRecurrenceCoeffs(const unsigned k) const
	{
	const Recurrence& rc(rCoeffs_.at(k));
	return std::pair<precise_type,precise_type>(rc.alpha, rc.beta);
	}

	private:
	static const precise_type precise_zero;

	struct Recurrence
	{
	inline Recurrence(const precise_type a,
	const precise_type b)
	: alpha(a), beta(b)
	{
	assert(beta > 0.0);
	sqrbeta = std::sqrt(beta);
	}

	precise_type alpha;
	precise_type beta;
	precise_type sqrbeta;
	};

	class PolyFcn;
	friend class PolyFcn;

	class PolyFcn : public Functor1<precise_type, precise_type>
	{
	public:
	inline PolyFcn(const ContOrthoPoly1D& p,
	const unsigned d1, const unsigned power)
	: poly(p), deg1(d1), polypow(power) {}

	inline virtual ~PolyFcn() {}

	inline precise_type operator()(const precise_type& x) const
	{
	precise_type p = 1.0;
	if (polypow)
	{
	p = poly.normpoly(deg1, x);
	switch (polypow)
	{
	case 1U:
	break;
	case 2U:
	p *= p;
	break;
	default:
	const precise_type myPow = polypow;
	p = std::pow(p, myPow);
	break;
	}
	}
	return p;
	}

	private:
	const ContOrthoPoly1D& poly;
	unsigned deg1;
	unsigned polypow;
	};

	class TripleProdFcn;
	friend class TripleProdFcn;

	class TripleProdFcn : public Functor1<precise_type, precise_type>
	{
	public:
	inline TripleProdFcn(const ContOrthoPoly1D& p, const unsigned d1,
	const unsigned d2, const unsigned d3)
	: poly(p)
	{
	degs[0] = d1;
	degs[1] = d2;
	degs[2] = d3;
	degmax = std::max(d1, std::max(d2, d3));
	}

	inline virtual ~TripleProdFcn() {}

	inline precise_type operator()(const precise_type& x) const
	{return poly.normpolyprod(degs, 3, degmax, x);}

	private:
	const ContOrthoPoly1D& poly;
	unsigned degs[3];
	unsigned degmax;
	};

	template<class Funct> class EWPolyProductFcn;
	template<class Funct> friend class EWPolyProductFcn;

	template<class Funct>
	class EWPolyProductFcn : public Functor1<precise_type, precise_type>
	{
	public:
	inline EWPolyProductFcn(const ContOrthoPoly1D& p, const Funct& w,
	const unsigned d1, const unsigned d2)
	: poly(p), wf(w)
	{
	degs[0] = d1;
	degs[1] = d2;
	degmax = std::max(d1, d2);
	}

	inline virtual ~EWPolyProductFcn() {}

	inline precise_type operator()(const precise_type& x) const
	{return poly.normpolyprod(degs,2,degmax,x-poly.shiftX_)*wf(x);}

	private:
	const ContOrthoPoly1D& poly;
	const Funct& wf;
	unsigned degs[2];
	unsigned degmax;
	};

	class MemoKey
	{
	public:
	inline MemoKey(const unsigned d0, const unsigned d1,
	const unsigned d2, const unsigned d3)
	: nDegs_(4U)
	{
	degs_[0] = d0;
	degs_[1] = d1;
	degs_[2] = d2;
	degs_[3] = d3;
	std::sort(degs_, degs_+nDegs_);
	}

	inline MemoKey(const unsigned d0, const unsigned d1,
	const unsigned d2, const unsigned d3,
	const unsigned d4, const unsigned d5)
	: nDegs_(6U)
	{
	degs_[0] = d0;
	degs_[1] = d1;
	degs_[2] = d2;
	degs_[3] = d3;
	degs_[4] = d4;
	degs_[5] = d5;
	std::sort(degs_, degs_+nDegs_);
	}

	inline MemoKey(const unsigned d0, const unsigned d1,
	const unsigned d2, const unsigned d3,
	const unsigned d4, const unsigned d5,
	const unsigned d6, const unsigned d7)
	: nDegs_(8U)
	{
	degs_[0] = d0;
	degs_[1] = d1;
	degs_[2] = d2;
	degs_[3] = d3;
	degs_[4] = d4;
	degs_[5] = d5;
	degs_[6] = d6;
	degs_[7] = d7;
	std::sort(degs_, degs_+nDegs_);
	}

	inline const unsigned* degrees() const {return degs_;}
	inline unsigned nDegrees() const {return nDegs_;}

	inline bool operator<(const MemoKey& r) const
	{
	if (nDegs_ < r.nDegs_)
	return true;
	if (r.nDegs_ < nDegs_)
	return false;
	for (unsigned i=0; i<nDegs_; ++i)
	{
	if (degs_[i] < r.degs_[i])
	return true;
	if (r.degs_[i] < degs_[i])
	return false;
	}
	return false;
	}

	inline bool operator==(const MemoKey& r) const
	{
	if (nDegs_ != r.nDegs_)
	return false;
	for (unsigned i=0; i<nDegs_; ++i)
	if (degs_[i] != r.degs_[i])
	return false;
	return true;
	}

	inline bool operator!=(const MemoKey& r) const
	{return !(*this == r);}

	private:
	unsigned degs_[8];
	unsigned nDegs_;
	};

	template<typename Numeric1, typename Numeric2>
	void copyMeasure(const std::vector<std::pair<Numeric1,Numeric2> >& inMeasure);

	template<typename Numeric>
	void copyInputCoords(const std::vector<Numeric>& inCoords);

	void calcMeanXandWeightSums(bool shiftTheMeasure);
	void calcRecurrenceCoeffs(OrthoPolyMethod m);

	void calcRecurrenceStieltjes();
	void calcRecurrenceLanczos();

	template<class OrthoPoly>
	long double sigma(const OrthoPoly& ops, int k, int l,
	Matrix<long double>& sigMatrix,
	Matrix<int>& flagMatrix) const;

	template<class OrthoPoly>
	void calculateMonicAverages(const OrthoPoly& ops, long double* averages,
	unsigned maxdeg);

	template<class OrthoPoly>
	void calcRecurrenceModifiedChebyshev(const OrthoPoly& ops);

	precise_type monicpoly(unsigned degree, precise_type x) const;

	std::pair<precise_type,precise_type> monicInnerProducts(
	unsigned degree) const;

	precise_type normpoly(unsigned degree, precise_type x) const;

	std::pair<precise_type,precise_type> twonormpoly(
	unsigned deg1, unsigned deg2, precise_type x) const;

	precise_type normpolyprod(const unsigned* degrees, unsigned nDeg,
	unsigned degmax, precise_type x) const;

	template <typename Numeric>
	precise_type normseries(const Numeric* coeffs, unsigned maxdeg,
	precise_type x) const;

	template <typename Numeric>
	void getAllPolys(double x, unsigned maxdeg, Numeric *polyValues) const;

	double getCachedAverage(const MemoKey& key) const;

	std::vector<MeasurePoint> measure_;
	std::vector<Recurrence> rCoeffs_;
	mutable std::map<MemoKey,double> cachedAverages_;
	precise_type wsum_;
	precise_type wsumsq_;
	precise_type meanX_;
	precise_type shiftX_;
	double minX_;
	double maxX_;
	unsigned maxdeg_;
	bool allWeightsEqual_;

	#ifdef SWIG
	public:
	inline StorablePolySeries1D* makeStorablePolySeries_2(
	double xmin, double xmax, const std::vector<double>& coeffs) const
	{
	CPP11_auto_ptr<StorablePolySeries1D> ptr;
	if (coeffs.empty())
	ptr = this->makeStorablePolySeries(xmin, xmax);
	else
	ptr = this->makeStorablePolySeries(xmin, xmax, &coeffs[0],
	coeffs.size() - 1);
	return ptr.release();
	}

	inline double sumOfWeights_2() const
	{return sumOfWeights();}

	inline double sumOfWeightSquares_2() const
	{return sumOfWeightSquares();}

	inline double meanCoordinate_2() const
	{return meanCoordinate();}

	inline double empiricalKroneckerDelta_2(const unsigned deg1,
	const unsigned deg2) const
	{return empiricalKroneckerDelta(deg1, deg2);}

	inline Matrix<double> jacobiMatrix_2(const unsigned n) const
	{return jacobiMatrix(n);}

	inline double integratePoly_2(const unsigned degree, const unsigned power,
	const double xmin, const double xmax) const
	{return integratePoly(degree, power, xmin, xmax);}

	inline double integrateTripleProduct_2(const unsigned deg1, const unsigned deg2,
	const unsigned deg3, const double xmin,
	const double xmax) const
	{return integrateTripleProduct(deg1, deg2, deg3, xmin, xmax);}

	inline Matrix<double> tripleProductMatrix_2(
	const std::vector<std::pair<unsigned,unsigned> >& pairs,
	const unsigned maxdeg, const double xmin, const double xmax) const
	{return tripleProductMatrix(pairs, maxdeg, xmin, xmax);}

	inline double location2() const {return location();}

	inline double scale2() const {return scale();}

	inline std::pair<double,double> monicRecurrenceCoeffs2(const unsigned deg) const
	{
	const std::pair<precise_type,precise_type>& p = monicRecurrenceCoeffs(deg);
	return std::pair<double,double>(p.first, p.second);
	}

	inline std::pair<double,double> recurrenceCoeffs2(const unsigned deg) const
	{
	const std::pair<precise_type,precise_type>& p = recurrenceCoeffs(deg);
	return std::pair<double,double>(p.first, p.second);
	}
	#endif // SWIG
	};

	/**
	// A functor for a particular degree polynomial of the given
	// ContOrthoPoly1D. The poly system is not copied, only a reference
	// is used. It is a responsibility of the user to make sure that the
	// lifetime of the poly system object exceeds the lifetime of the functor.
	*/
	class ContOrthoPoly1DDeg : public Functor1<long double, long double>
	{
	public:
	inline ContOrthoPoly1DDeg(const ContOrthoPoly1D& fcn,
	const unsigned degree,
	const long double normfactor=1.0L)
	: fcn_(fcn), norm_(normfactor), deg_(degree)
	{
	if (deg_ > fcn_.maxDegree()) throw std::invalid_argument(
	"In npstat::ContOrthoPoly1DDeg::constructor: "
	"degree argument is out of range");
	}

	inline virtual ~ContOrthoPoly1DDeg() {}

	inline virtual long double operator()(const long double& a) const
	{return norm_*fcn_.poly(deg_, a);}

	private:
	ContOrthoPoly1DDeg();

	const ContOrthoPoly1D& fcn_;
	long double norm_;
	unsigned deg_;
	};
	}

	#include "npstat/nm/ContOrthoPoly1D.icc"

	#endif // NPSTAT_CONTORTHOPOLY1D_HH_

File Metadata

Mime Type: text/x-c
Expires: Tue, Sep 30, 4:45 AM (15 h, 38 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 6564809
Default Alt Text: ContOrthoPoly1D.hh (31 KB)

ContOrthoPoly1D.hhNo OneTemporaryActions

ContOrthoPoly1D.hhView Options

File Metadata

Event Timeline

ContOrthoPoly1D.hh
No OneTemporary
Actions

ContOrthoPoly1D.hh
View Options