KernelSensitivityCalculator.hh

Go to the documentation of this file.

#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 <utility>
 
#include "npstat/nm/Matrix.hh"
#include "npstat/nm/AbsIntervalQuadrature1D.hh"
 
#include "npstat/stat/AbsDistribution1D.hh"
 
namespace npstat {
    class KernelSensitivityCalculator
    {
    public:
        /**
        // 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 <class Kernel>
        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, bool validateCDF = true);
 
        /**
        // 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 <class Kernel>
        KernelSensitivityCalculator(
            const Kernel& kernel, const std::vector<double>& sample,
            unsigned maxdegInputPoly, unsigned maxdegOutPoly,
            long double ymin, long double ymax,
            const AbsIntervalQuadrature1D& yIntegrator,
            bool normalizeKernel = true, bool validateCDF = true);
 
        inline bool isOracle() const {return oracle_;}
        inline unsigned maxdegInputPoly() const {return maxdegInputPoly_;}
        inline unsigned maxdegOutPoly() const {return maxdegOutPoly_;}
 
        /**
        // Return the "standard" kernel sensitivity matrix.
        //
        // "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 orthonormal with the weight
        // given by the KDE density.
        */
        template<class InPoly, class OutPoly>
        Matrix<double> sMatrix0(
            const InPoly& inputPoly, const OutPoly& outPoly) const;
 
        /**
        // Return the kernel sensitivity matrix that uses
        // comparison density in y.
        //
        // "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<class InPoly, class KDECdf>
        Matrix<double> sMatrix1(
            const InPoly& inputPoly, const KDECdf& cdf) const;
 
        /**
        // Return the kernel sensitivity matrix that uses
        // comparison density in x.
        //
        // "outPoly" is the OPS orthonormal with the weight
        // given by the KDE density.
        */
        template<class OutPoly>
        Matrix<double> 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<class KDECdf>
        Matrix<double> sMatrix3(const KDECdf& cdf) const;
 
    private:
        typedef std::vector<std::pair<long double,long double> > WeightedCoords;
 
        template <class Kernel>
        void calcKValues(const Kernel& kernel, bool normalize);
 
        void makeLocalLegendreMatrix() const;
 
        template<class Poly>
        Matrix<long double> standardPolyMatrix(
            const Poly& poly, unsigned maxdeg, const WeightedCoords& wcoords) const;
 
        template<class CdfFcn>
        Matrix<long double> legendrePolyMatrix(
            const CdfFcn& cdf, unsigned maxdeg, const WeightedCoords& wcoords) const;
 
        std::shared_ptr<const AbsDistribution1D> distro_;
        unsigned maxdegInputPoly_;
        unsigned maxdegOutPoly_;
        Matrix<long double> kValues_;
        WeightedCoords xw_;
        WeightedCoords yw_;
        mutable Matrix<long double> localLegendreMatrix_;
        long double densityIntegral_;
        bool oracle_;
        bool validateCdf_;
    };
}
 
#include "npstat/stat/KernelSensitivityCalculator.icc"
 
#endif // NPSTAT_KERNELSENSITIVITYCALCULATOR_HH_