KDEFilterND.hh

Go to the documentation of this file.

#ifndef NPSTAT_KDEFILTERND_HH_
#define NPSTAT_KDEFILTERND_HH_
 
/*!
// \file KDEFilterND.hh
//
// \brief KDE on a regularly spaced multivariate grid with symmetric kernels
//
// If you don't expect to perform cross-validation or bandwidth scans,
// ConstantBandwidthSmootherND class might be a more convenient choice.
//
// Author: I. Volobouev
//
// December 2011
*/
 
#include <algorithm>
 
#include "npstat/stat/AbsPolyFilterND.hh"
 
#include "npstat/nm/ArrayND.hh"
#include "npstat/nm/ConvolutionEngineND.hh"
 
namespace npstat {
    /**
    // This class performs kernel density estimation on a regularly spaced
    // multivariate grid with arbitrary, tabulated symmetric kernels
    */
    template <unsigned MaxDeg>
    class KDEFilterND : public AbsPolyFilterND
    {
    public:
        enum {deg_size = MaxDeg};
 
        /**
        // Main constructor. The arguments are as follows:
        //
        //  taper      -- Damping factors for each polynomial degree
        //                (starting with the 0th order term). This can be
        //                NULL in which case it is assumed that all
        //                factors are 1.
        //
        //  maxDegree  -- Maximum degree of the polynomials. The value
        //                must not exceed the "MaxDeg" template parameter.
        //                The length of the "taper" array (if not NULL)
        //                must be equal to maxDegree + 1. Note that the same
        //                filter will be produced using an even degree N and
        //                an odd degree N+1.
        //
        //  weight     -- The array of weights. By default, this is just 
        //                one hyperoctant corresponding to positive directions
        //                in all dimensions. The remaining part of the kernel
        //                is generated by reflections. If the full weight is
        //                provided instead, change the "mirrorWeight"
        //                default.
        //
        //  engine     -- The engine to use for convolutions. This object
        //                will not make its own copy of the engine. It is
        //                the responsibility of the user to ensure that
        //                the lifetime of the engine exceeds the lifetime
        //                of this object. It is assumed that the size
        //                of the arrays to filter will be two times smaller
        //                in each dimension than what is expected by the
        //                engine. Note that dimensionalities of weight
        //                array must not exceed data dimensionalities.
        //
        //  id         -- Filter slot number in the engine. Each filter
        //                used with this particular engine must be provided
        //                with a unique slot id.
        //
        //  workBuffer -- Buffer array for various calculations. Must have
        //                the same shape as the convolution engine. Will
        //                not be copied internally. Must exist while this
        //                object operates.
        //
        //  mirrorData -- If true, the data will be mirrored at the
        //                boundaries to reduce the edge effects.
        //
        //  mirrorWeight -- If true, the weight will be mirrored. If false,
        //                  the weight is considered to be complete. In the
        //                  latter case the array size in each dimension
        //                  must be odd and the weight should be centered
        //                  in the array.
        //
        //  manageSlot -- If true, the Fourier transform of this filter
        //                will be removed from the engine filter collection
        //                in the destructor of this object.
        */
        template <class Array>
        KDEFilterND(const double* taper, unsigned maxDegree,
                    const Array& weight, ConvolutionEngineND& engine,
                    unsigned long id, ArrayND<double>& workBuffer,
                    bool mirrorData=true, bool mirrorWeight=true,
                    bool manageSlot=true);
 
        virtual ~KDEFilterND();
 
        //@{
        /** Inspect object properties */
        inline unsigned dim() const {return dataShape_.size();}
        inline unsigned maxDegree() const {return taper_.size() - 1U;}
        inline ArrayShape dataShape() const {return dataShape_;}
        double taper(unsigned degree) const;
        inline bool mirrorsData() const {return mirrorData_;}
        inline unsigned long engineSlotId() const {return slotId_;}
        inline double filterAtTheCenter() const {return filterAtTheCenter_;}
        //@}
 
        //@{
        /**
        // Contribution of a single point into the density estimate
        // at that point (not normalized). This is needed for various
        // leaving-one-out cross validation procedures.
        */
        double selfContribution(const unsigned* index,
                                unsigned indexLen) const;
        double linearSelfContribution(unsigned long index) const;
        //@}
 
        /** Check if the self contribution varies from point to point */
        inline bool selfContributionVaries() const {return selfContrib_;}
 
        /** Check the compatibility of an array with the filter */
        template <class Array>
        inline bool isCompatible(const Array& in) const
            {return in.isCompatible(dataShape_);}
 
        /** The method which performs the filtering */
        template <class ArrIn, class ArrOut>
        void filter(const ArrIn& in, ArrOut* out) const;
 
        /**
        // This method is included for compatibility with LocalPolyFilterND.
        // For KDE, there is no difference between "convolve" and "filter".
        */
        template <class ArrIn, class ArrOut>
        inline void convolve(const ArrIn& in, ArrOut* out) const
            {this->filter(in, out);}
 
        static inline unsigned classMaxDegree() {return MaxDeg;}
 
    private:
        KDEFilterND();
        KDEFilterND(const KDEFilterND&);
        KDEFilterND& operator=(const KDEFilterND&);
 
        template <class Array>
        void initialize(const double* itaper, unsigned maxDegree,
                        const Array& inweight, bool mirrorWeight);
 
        void calculateSelfContribution(const ArrayND<double>& filter);
 
        ArrayShape dataShape_;
        mutable ArrayShape idxBuffer_;
        ConvolutionEngineND* engine_;
        std::vector<double> taper_;
        unsigned long slotId_;
        ArrayND<double>* buf_;
        double filterAtTheCenter_;
        ArrayND<double>* selfContrib_;
        bool mirrorData_;
        bool manageSlot_;
 
#ifdef SWIG
    public:
        inline KDEFilterND(const std::vector<double>& itaper,
                           const ArrayND<double>& inweight,
                           ConvolutionEngineND& engine,
                           unsigned long id, ArrayND<double>& workBuffer,
                           bool mirrorData=true, bool mirrorWeight=true,
                           bool manageSlot=true)
            : dataShape_(halfShape(engine.shape())),
              idxBuffer_(inweight.shape()),
              engine_(&engine),
              taper_(std::max(itaper.size(), static_cast<std::size_t>(1)), 1.0),
              slotId_(id),
              buf_(&workBuffer),
              filterAtTheCenter_(0.0),
              selfContrib_(0),
              mirrorData_(mirrorData),
              manageSlot_(manageSlot)
        {
            const double* tap = 0;
            if (!itaper.empty())
                tap = &itaper[0];
            initialize(tap, taper_.size()-1, inweight, mirrorWeight);
        }
#endif //SWIG
    };
}
 
#include "npstat/stat/KDEFilterND.icc"
 
#endif // NPSTAT_KDEFILTERND_HH_