Index: trunk/npstat/stat/KernelSensitivityCalculator.hh =================================================================== --- trunk/npstat/stat/KernelSensitivityCalculator.hh (revision 892) +++ trunk/npstat/stat/KernelSensitivityCalculator.hh (revision 893) @@ -1,124 +1,179 @@ #ifndef NPSTAT_KERNELSENSITIVITYCALCULATOR_HH_ #define NPSTAT_KERNELSENSITIVITYCALCULATOR_HH_ /*! // \file KernelSensitivityCalculator.hh // // \brief Helper class for calculating kernel sensitivity matrices // // Author: I. Volobouev // // March 2023 */ #include #include "npstat/nm/Matrix.hh" #include "npstat/nm/AbsIntervalQuadrature1D.hh" #include "npstat/stat/AbsDistribution1D.hh" namespace npstat { class KernelSensitivityCalculator { public: - /** Oracle constructor */ + /** + // Oracle constructor. Arguments are as follows: + // + // kernel The kernel functor, K(x,y). Must have a method + // with a signature similar to + // "double operator()(double x, double y) const". + // + // distro The oracle distribution. + // + // maxdegInputPoly Maximum degree of P_k(x) polynomials. + // + // maxdegOutPoly Maximum degree of Q_j(y) polynomials. + // + // ymin, ymax The interval on which the KDE of the oracle + // density is supported. + // + // xIntegrator Numerical quadrature object for calculating + // integrals in x. + // + // yIntegrator Numerical quadrature object for calculating + // integrals in y. + // + // normalizeKernel Set this to "true" in order to normalize the + // kernel internally so that Integral K(x,y) dy = 1 + // for every x. Set this to "false" if the kernel + // is already normalized. + */ template KernelSensitivityCalculator( const Kernel& kernel, const AbsDistribution1D& distro, unsigned maxdegInputPoly, unsigned maxdegOutPoly, long double ymin, long double ymax, const AbsIntervalQuadrature1D& xIntegrator, const AbsIntervalQuadrature1D& yIntegrator, bool normalizeKernel = true); - /** Non-oracle constructor */ + /** + // Non-oracle constructor. Arguments are as follows: + // + // kernel The kernel functor, K(x,y). Must have a method + // with a signature similar to + // "double operator()(double x, double y) const". + // + // sample The sample of points on which to perform KDE. + // + // maxdegInputPoly Maximum degree of P_k(x) polynomials. + // + // maxdegOutPoly Maximum degree of Q_j(y) polynomials. + // + // ymin, ymax The interval on which the KDE is supported. + // + // yIntegrator Numerical quadrature object for calculating + // integrals in y. + // + // normalizeKernel Set this to "true" in order to normalize the + // kernel internally so that Integral K(x,y) dy = 1 + // for every x. Set this to "false" if the kernel + // is already normalized. + // + // validateCDF If "true", the code will check that all values + // of the KDE cumulative distribution function + // belong to the [0, 1] interval and will throw + // an exception if this is not the case. If "false", + // this check will not be performed. + */ template KernelSensitivityCalculator( const Kernel& kernel, const std::vector& sample, unsigned maxdegInputPoly, unsigned maxdegOutPoly, long double ymin, long double ymax, const AbsIntervalQuadrature1D& yIntegrator, - bool normalizeKernel = true); + bool normalizeKernel = true, bool validateCDF = true); inline bool isOracle() const {return oracle_;} /** // Return the "standard" kernel sensitivity matrix. // - // "inputPoly" is the OPS orthogonal with the weight + // "inputPoly" is the OPS orthonormal with the weight // "distro" (oracle case) or with the discrete measure // generated by the sample (non-oracle case). // - // "outPoly" is the OPS orthogonal with the weight + // "outPoly" is the OPS orthonormal with the weight // given by the KDE density. */ template Matrix sMatrix0( const InPoly& inputPoly, const OutPoly& outPoly) const; /** // Return the kernel sensitivity matrix that uses // comparison density in y. // - // "inputPoly" is the OPS orthogonal with the weight + // "inputPoly" is the OPS orthonormal with the weight // "distro" (oracle case) or with the discrete measure // generated by the sample (non-oracle case). // // "cdf" is the functor that returns the cumulative // distribution function for the KDE density. */ template Matrix sMatrix1( const InPoly& inputPoly, const KDECdf& cdf) const; /** // Return the kernel sensitivity matrix that uses // comparison density in x. // - // "outPoly" is the OPS orthogonal with the weight + // "outPoly" is the OPS orthonormal with the weight // given by the KDE density. */ template Matrix sMatrix2(const OutPoly& outPoly) const; /** // Return the kernel sensitivity matrix that uses // comparison density in both x and y. // // "cdf" is the functor that returns the cumulative // distribution function for the KDE density. */ template Matrix sMatrix3(const KDECdf& cdf) const; private: typedef std::vector > WeightedCoords; template void calcKValues(const Kernel& kernel, bool normalize); void makeLocalLegendreMatrix() const; template Matrix standardPolyMatrix( const Poly& poly, unsigned maxdeg, const WeightedCoords& wcoords) const; template Matrix legendrePolyMatrix( const CdfFcn& cdf, unsigned maxdeg, const WeightedCoords& wcoords) const; std::shared_ptr distro_; unsigned maxdegInputPoly_; unsigned maxdegOutPoly_; Matrix kValues_; WeightedCoords xw_; WeightedCoords yw_; mutable Matrix localLegendreMatrix_; bool oracle_; + bool validateCdf_; }; } #include "npstat/stat/KernelSensitivityCalculator.icc" #endif // NPSTAT_KERNELSENSITIVITYCALCULATOR_HH_ Index: trunk/npstat/stat/KernelSensitivityCalculator.icc =================================================================== --- trunk/npstat/stat/KernelSensitivityCalculator.icc (revision 892) +++ trunk/npstat/stat/KernelSensitivityCalculator.icc (revision 893) @@ -1,201 +1,203 @@ #include #include #include #include #include #include #include "npstat/nm/SimpleFunctors.hh" #include "npstat/nm/ClassicalOrthoPolys1D.hh" namespace npstat { template KernelSensitivityCalculator::KernelSensitivityCalculator( const Kernel& kernel, const AbsDistribution1D& distro, const unsigned maxdegInputPoly, const unsigned maxdegOutPoly, const long double ymin, const long double ymax, const AbsIntervalQuadrature1D& xIntegrator, const AbsIntervalQuadrature1D& yIntegrator, const bool normalizeKernel) : distro_(distro.clone()), maxdegInputPoly_(maxdegInputPoly), maxdegOutPoly_(maxdegOutPoly), kValues_(yIntegrator.npoints(), xIntegrator.npoints()), xw_(xIntegrator.weightedIntegrationPointsLD( DensityFunctor1D(distro), distro.quantile(0.0), distro.quantile(1.0))), yw_(yIntegrator.weightedIntegrationPointsLD( ConstValue1(1.0L), ymin, ymax)), - oracle_(true) + oracle_(true), + validateCdf_(true) { // Renormalize the density long double densityIntegral = 0.0L; const unsigned long sz = xw_.size(); for (unsigned long i=0; i KernelSensitivityCalculator::KernelSensitivityCalculator( const Kernel& kernel, const std::vector& sample, const unsigned maxdegInputPoly, const unsigned maxdegOutPoly, const long double ymin, const long double ymax, const AbsIntervalQuadrature1D& yIntegrator, - const bool normalizeKernel) + const bool normalizeKernel, const bool validateCDF) : maxdegInputPoly_(maxdegInputPoly), maxdegOutPoly_(maxdegOutPoly), kValues_(yIntegrator.npoints(), sample.size()), xw_(sample.size()), yw_(yIntegrator.weightedIntegrationPointsLD( ConstValue1(1.0L), ymin, ymax)), - oracle_(false) + oracle_(false), + validateCdf_(validateCDF) { const unsigned long sz = sample.size(); if (sz == 0UL || sz < maxdegInputPoly*2U) throw std::invalid_argument( "In npstat::KernelSensitivityCalculator constructor: " "insufficient sample size"); const long double w = 1.0L/sz; for (unsigned long i=0; i(sample[i], w); std::sort(xw_.begin(), xw_.end()); // Calculate all necessary kernel values calcKValues(kernel, normalizeKernel); } template void KernelSensitivityCalculator::calcKValues( const Kernel& kernel, const bool normalize) { const unsigned nx = kValues_.nColumns(); assert(nx == xw_.size()); const unsigned ny = kValues_.nRows(); assert(ny == yw_.size()); assert(maxdegInputPoly_); assert(maxdegOutPoly_); if (normalize) { for (unsigned ix=0; ix 0.0L); for (unsigned iy=0; iy Matrix KernelSensitivityCalculator::standardPolyMatrix( const Poly& poly, const unsigned maxdeg, const WeightedCoords& wcoords) const { const unsigned nPoints = wcoords.size(); Matrix mat(nPoints, maxdeg); std::vector buf(maxdeg + 1U); for (unsigned ipt=0; ipt Matrix KernelSensitivityCalculator::legendrePolyMatrix( const CdfFcn& cdf, const unsigned maxdeg, const WeightedCoords& wcoords) const { const unsigned nPoints = wcoords.size(); Matrix mat(nPoints, maxdeg); std::vector buf(maxdeg + 1U); ShiftedLegendreOrthoPoly1D poly; for (unsigned ipt=0; ipt 1.0L) + if (validateCdf_ && (cdfValue < 0.0L || cdfValue > 1.0L)) { std::ostringstream os; - os.precision(std::numeric_limits::digits10); + os.precision(std::numeric_limits::digits10 + 1); os << "In npstat::KernelSensitivityCalculator::legendrePolyMatrix: " << "cdf value " << cdfValue << " at " << wcoords[ipt].first << " is out of range"; throw std::runtime_error(os.str()); } poly.allpoly(cdfValue, &buf[0], maxdeg); long double* row = mat.data() + maxdeg*ipt; const long double* from = &buf[1]; for (unsigned deg=1; deg<=maxdeg; ++deg) *row++ = *from++; } return mat; } template Matrix KernelSensitivityCalculator::sMatrix0( const InPoly& inputPoly, const OutPoly& outPoly) const { const Matrix& pMatrix = standardPolyMatrix( inputPoly, maxdegInputPoly_, xw_); const Matrix& qMatrix = standardPolyMatrix( outPoly, maxdegOutPoly_, yw_); return qMatrix.T()*kValues_*pMatrix; } template Matrix KernelSensitivityCalculator::sMatrix1( const InPoly& inputPoly, const KDECdf& cdf) const { const Matrix& pMatrix = standardPolyMatrix( inputPoly, maxdegInputPoly_, xw_); const Matrix& qMatrix = legendrePolyMatrix( cdf, maxdegOutPoly_, yw_); return qMatrix.T()*kValues_*pMatrix; } template Matrix KernelSensitivityCalculator::sMatrix2( const OutPoly& outPoly) const { makeLocalLegendreMatrix(); const Matrix& qMatrix = standardPolyMatrix( outPoly, maxdegOutPoly_, yw_); return qMatrix.T()*kValues_*localLegendreMatrix_; } template Matrix KernelSensitivityCalculator::sMatrix3( const KDECdf& cdf) const { makeLocalLegendreMatrix(); const Matrix& qMatrix = legendrePolyMatrix( cdf, maxdegOutPoly_, yw_); return qMatrix.T()*kValues_*localLegendreMatrix_; } }