LocalPolyFilterND.hh

Go to the documentation of this file.

#ifndef NPSTAT_LOCALPOLYFILTERND_HH_
#define NPSTAT_LOCALPOLYFILTERND_HH_
 
/*!
// \file LocalPolyFilterND.hh
//
// \brief Local polynomial filtering (regression) on uniform
//        hyperrectangular grids
//
// Author: I. Volobouev
//
// November 2009
*/
 
#include <map>
#include <vector>
 
#include "geners/ClassId.hh"
#include "geners/CPP11_auto_ptr.hh"
 
#include "npstat/nm/ArrayND.hh"
#include "npstat/nm/Matrix.hh"
 
#include "npstat/stat/AbsPolyFilterND.hh"
 
namespace npstat {
    typedef ArrayND<double> PolyFilterND;
 
    /** This class performs local polynomial filtering in multiple dimensions */
    template <unsigned MaxDeg>
    class LocalPolyFilterND : 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, far
        //               away from the boundaries (where the situation is
        //               symmetric) the same filter will be produced using
        //               the same taper with an even degree N and with
        //               an odd degree N+1. Near the boundaries the filter
        //               coefficients will, of course, differ in these
        //               two cases.
        //
        //  weight    -- The array of weights.
        //
        //  mirrorWeight -- If true, the "weight" argument represents
        //                  just one hyperoctant corresponding to
        //                  positive directions in all dimensions.
        //                  The remaining part of the weight should
        //                  be generated by reflections. If false,
        //                  the "weight" argument will not be mirrored.
        //                  In the latter case, the maximum of the weight
        //                  should be in the array center.
        //
        //  dataShape -- The shape of the data arrays which will be
        //               processed by this filter
        */
        template <class Array>
        LocalPolyFilterND(const double* taper, unsigned maxDegree,
                          const Array& weight, bool mirrorWeight,
                          const ArrayShape& dataShape);
 
        LocalPolyFilterND(const LocalPolyFilterND&);
        LocalPolyFilterND& operator=(const LocalPolyFilterND&);
 
        virtual ~LocalPolyFilterND();
 
        //@{
        /** Inspect object properties */
        inline unsigned dim() const {return bins_.rank();}
        inline unsigned maxDegree() const {return maxDegree_;}
        inline ArrayShape dataShape() const {return bins_.shape();}
        double taper(unsigned degree) const;
        //@}
 
        /** Get the filter coefficients for the given grid point */
        const PolyFilterND& getFilter(const unsigned* index,
                                      unsigned lenIndex) const;
 
        /** Extract filter coefficients using linear grid index */
        const PolyFilterND& linearGetFilter(unsigned long index) const;
 
        /** Get the complete effective filter matrix */
        Matrix<double> getFilterMatrix() const;
 
        /** Get the info needed to construct the sparse filter matrix */
        template <class Triplet>
        CPP11_auto_ptr<std::vector<Triplet> > sparseFilterTriplets() const;
 
        //@{
        /**
        // 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 lenIndex) const;
        double linearSelfContribution(unsigned long index) const;
        //@}
 
        /** Check compatibility of an array with the filter */
        template <class Array>
        inline bool isCompatible(const Array& in) const
            {return bins_.isShapeCompatible(in);}
 
        /** This method performs the filtering */
        template <class ArrIn, class ArrOut>
        void filter(const ArrIn& in, ArrOut* out) const;
 
        /**
        // A diffent filtering method in which the shapes of the
        // kernels are determined by the positions of the "sources"
        // (i.e., sample points) instead of the positions at which
        // the density (or response) is estimated. Note that elements
        // of "out" array themselves are used as result accumulators.
        */
        template <class ArrIn, class ArrOut>
        void convolve(const ArrIn& in, ArrOut* out) const;
 
        template <unsigned MaxDeg2>
        bool operator==(const LocalPolyFilterND<MaxDeg2>& r) const;
 
        template <unsigned MaxDeg2>
        inline bool operator!=(const LocalPolyFilterND<MaxDeg2>& r) const
            {return !(*this == r);}
 
        // Methods needed for I/O
        inline virtual gs::ClassId classId() const {return gs::ClassId(*this);}
        virtual bool write(std::ostream& os) const;
 
        static const char* classname();
        static inline unsigned version() {return 1;}
        static LocalPolyFilterND* read(const gs::ClassId& id, std::istream& in);
 
        static inline unsigned classMaxDegree() {return MaxDeg;}
 
    private:
        typedef std::map<ArrayRange, PolyFilterND*> PolyMap;
 
        LocalPolyFilterND();
 
        template <class Array>
        void initialize(const double* itaper, unsigned maxDegree,
                        const Array& inweight, bool mirrorWeight);
 
        void cleanup();
        void copyFilters(const LocalPolyFilterND& r);
 
        double* taper_;
        PolyMap unique_;
        ArrayShape wshape_;
        ArrayND<PolyFilterND*> bins_;
        std::vector<unsigned> filterCenters_;
        PolyFilterND* centerFilter_;
        unsigned maxDegree_;
        mutable std::vector<unsigned> indexBuf_;
 
        template <unsigned StackLen, unsigned StackDim>
        PolyFilterND* buildFilter(const ArrayND<double,StackLen,StackDim>& w,
                                  const unsigned* center,
                                  const double* steps, unsigned nSteps);
 
        template <unsigned StackLen, unsigned StackDim>
        void createFiltersLoop(unsigned level, unsigned long idxData,
                               const ArrayND<double,StackLen,StackDim>& weight,
                               ArrayRange& range, unsigned* center,
                               const double* steps, unsigned nSteps);
 
        template <typename Tin, unsigned StackLen, unsigned StackDim,
                  typename Tout, unsigned StackLen2, unsigned StackDim2>
        void filterLoop(unsigned level, unsigned long idxData,
                        const ArrayND<Tin,StackLen,StackDim>& in,
                        ArrayND<Tout,StackLen2,StackDim2>* out,
                        unsigned* shift) const;
 
        template <typename Tin, unsigned StackLen, unsigned StackDim>
        long double filterInnerLoop(unsigned level,
                                    const ArrayND<Tin,StackLen,StackDim>& in,
                                    const PolyFilterND& poly,
                                    const unsigned* shift,
                                    unsigned long idxData,
                                    unsigned long idxPoly) const;
 
        void matrixLoop(unsigned level, unsigned long idxData,
                        Matrix<double>* out, unsigned* shift) const;
 
        void matrixInnerLoop(unsigned level, Matrix<double>* out,
                             unsigned long rowNumber,
                             const PolyFilterND& poly, const unsigned* shift,
                             unsigned long idxData, unsigned long idxPoly) const;
 
        template <class Triplet>
        void sparseMatrixLoop(unsigned level, unsigned long idxData,
                              std::vector<Triplet>* out, unsigned* shift) const;
 
        template <class Triplet>
        void sparseMatrixInnerLoop(unsigned level, std::vector<Triplet>* out,
                             unsigned long rowNumber,
                             const PolyFilterND& poly, const unsigned* shift,
                             unsigned long idxData, unsigned long idxPoly) const;
 
        template <class ArrOut>
        void runConvolution(double w, unsigned long linearBin,
                            ArrOut* out, unsigned* cornerBuf, unsigned* indexBuf,
                            const unsigned* zeroBuf) const;
 
        static std::string generateClassName();
 
#ifdef SWIG
    public:
        inline LocalPolyFilterND(const double* itaper, unsigned maxDegree,
                                 const ArrayND<double>& inweight,
                                 bool mirrorWeight,
                                 const ArrayShape& dataShape)
            : taper_(0), wshape_(inweight.shape()), bins_(dataShape),
              centerFilter_(0), maxDegree_(maxDegree)
        {
            initialize(itaper, maxDegree, inweight, mirrorWeight);
        }
#endif //SWIG
    };
}
 
#include "npstat/stat/LocalPolyFilterND.icc"
 
#endif // NPSTAT_LOCALPOLYFILTERND_HH_