AbsBandwidthGCV.hh

Go to the documentation of this file.

#ifndef NPSTAT_ABSBANDWIDTHGCV_HH_
#define NPSTAT_ABSBANDWIDTHGCV_HH_
 
/*!
// \file AbsBandwidthGCV.hh 
//
// \brief Interface definitions for cross-validation with grouped data
//
// Implementations of these interfaces should calculate the quantity
// which is to be maximized (pseudo log likelihood, -AMISE, etc.)
//
// Author: I. Volobouev
//
// June 2015
*/
 
#include "npstat/stat/HistoND.hh"
 
#ifdef SWIG
#include <stdexcept>
#include "npstat/wrap/arrayNDToNumpy.hh"
#include "npstat/wrap/numpyArrayUtils.hh"
#endif // SWIG
 
namespace npstat {
    // Forward declarations
    struct AbsPolyFilter1D;
    struct AbsPolyFilterND;
 
    /**
    // Cross-validation for univariate densities.
    // "Numeric" template parameter is the type of histogram bins.
    // "Num2" is the type produced by the density estimation code
    // (usually just double).
    */
    template<typename Numeric, typename Num2>
    struct AbsBandwidthGCV1D
    {
        typedef Numeric bin_type;
        typedef Num2 density_type;
 
        inline virtual ~AbsBandwidthGCV1D() {}
 
        /**
        // It should be assumed that the "nFillsInRange" method of the
        // argument histogram returns the actual number of fills (that is,
        // the histogram represents an actual collection of points, has
        // possible bin values of 0, 1, 2, ..., and it is not scaled).
        */
        virtual double operator()(
            const HistoND<Numeric>& histo,
            const Num2* densityEstimate,
            const Num2* leaveOneOutEstimate,
            unsigned lenEstimate,
            const AbsPolyFilter1D& filterUsed) const = 0;
 
        /** Cross-validation for samples of univariate weighted points */
        virtual double operator()(
            const HistoND<Numeric>& histo,
            double effectiveSampleSize,
            const Num2* densityEstimate,
            const Num2* leaveOneOutEstimate,
            unsigned lenEstimate,
            const AbsPolyFilter1D& filterUsed) const = 0;
 
#ifdef SWIG
        inline double cv(
            const HistoND<Numeric,HistoAxis>& histo,
            PyObject* input, PyObject* leaveOneOut,
            const AbsPolyFilter1D& filterUsed) const
        {
            const int typenum = numpyArrayType(input);
            const int typenum2 = numpyArrayType(leaveOneOut);
            if (typenum != NumpyTypecode<double>::code ||
                typenum2 != NumpyTypecode<double>::code)
                throw std::invalid_argument("In npstat::AbsBandwidthGCV1D::cv: "
                                            "unsupported input array type");
            const unsigned datalen = histo.nBins();
            if (!isNumpyShapeCompatible(input, &datalen, 1U) ||
                !isNumpyShapeCompatible(leaveOneOut, &datalen, 1U)) throw std::invalid_argument(
                "In npstat::AbsBandwidthCV1D::cv: incompatible array shape");
            const Num2* in = (double*)PyArray_DATA((PyArrayObject*)input);
            const Num2* leave = (double*)PyArray_DATA((PyArrayObject*)leaveOneOut);
            return this->operator()(histo, in, leave, datalen, filterUsed);
        }
 
        inline double cvWeighted(
            const HistoND<Numeric,HistoAxis>& histo,
            const double effectiveSampleSize,
            PyObject* input, PyObject* leaveOneOut,
            const AbsPolyFilter1D& filterUsed) const
        {
            const int typenum = numpyArrayType(input);
            const int typenum2 = numpyArrayType(leaveOneOut);
            if (typenum != NumpyTypecode<double>::code ||
                typenum2 != NumpyTypecode<double>::code)
                throw std::invalid_argument("In npstat::AbsBandwidthGCV1D::cvWeighted: "
                                            "unsupported input array type");
            const unsigned datalen = histo.nBins();
            if (!isNumpyShapeCompatible(input, &datalen, 1U) ||
                !isNumpyShapeCompatible(leaveOneOut, &datalen, 1U)) throw std::invalid_argument(
                "In npstat::AbsBandwidthGCV1D::cvWeighted: incompatible array shape");
            const Num2* in = (double*)PyArray_DATA((PyArrayObject*)input);
            const Num2* leave = (double*)PyArray_DATA((PyArrayObject*)leaveOneOut);
            return this->operator()(histo, effectiveSampleSize, in,
                                    leave, datalen, filterUsed);
        }
#endif // SWIG
    };
 
    /**
    // Cross-validation for multivariate densities.
    // "Numeric" template parameter is the type of histogram bins.
    // "Array" template parameter should be one of ArrayND types.
    */
    template<typename Numeric, class Array>
    struct AbsBandwidthGCVND
    {
        typedef Numeric bin_type;
        typedef Array density_type;
 
        inline virtual ~AbsBandwidthGCVND() {}
 
        /**
        // It should be assumed that the "nFillsInRange" method of the
        // argument histogram returns the actual number of fills.
        */
        virtual double operator()(
            const HistoND<Numeric>& histo,
            const Array& densityEstimate,
            const Array& leaveOneOutEstimate,
            const AbsPolyFilterND& filterUsed) const = 0;
 
        /** Cross-validation for samples of multivariate weighted points */
        virtual double operator()(
            const HistoND<Numeric>& histo,
            double effectiveSampleSize,
            const Array& densityEstimate,
            const Array& leaveOneOutEstimate,
            const AbsPolyFilterND& filterUsed) const = 0;
 
#ifdef SWIG
        inline double cv(
            const HistoND<Numeric,HistoAxis>& histo,
            PyObject* input, PyObject* leaveOneOut,
            const AbsPolyFilterND& filterUsed) const
        {
            const int typenum = numpyArrayType(input);
            const int typenum2 = numpyArrayType(leaveOneOut);
            if (typenum != NumpyTypecode<double>::code ||
                typenum2 != NumpyTypecode<double>::code)
                throw std::invalid_argument("In npstat::AbsBandwidthCVGND::cv: "
                                            "unsupported input array type");
            const unsigned* shape = histo.binContents().shapeData();
            if (!isNumpyShapeCompatible(input, shape, histo.dim()) ||
                !isNumpyShapeCompatible(leaveOneOut, shape, histo.dim()))
                throw std::invalid_argument("In npstat::AbsBandwidthCVGND::cv: "
                                            "incompatible array shape");
            double* in = (double*)PyArray_DATA((PyArrayObject*)input);
            double* leave = (double*)PyArray_DATA((PyArrayObject*)leaveOneOut);
            const ArrayND<double>& wrap = externalMemArrayND(in, shape, histo.dim());
            const ArrayND<double>& wrapL = externalMemArrayND(leave, shape, histo.dim());
            return this->operator()(histo, wrap, wrapL, filterUsed);
        }
 
        inline double cvWeighted(
            const HistoND<Numeric,HistoAxis>& histo,
            const double effectiveSampleSize,
            PyObject* input, PyObject* leaveOneOut,
            const AbsPolyFilterND& filterUsed) const
        {
            const int typenum = numpyArrayType(input);
            const int typenum2 = numpyArrayType(leaveOneOut);
            if (typenum != NumpyTypecode<double>::code ||
                typenum2 != NumpyTypecode<double>::code)
                throw std::invalid_argument("In npstat::AbsBandwidthCVGND::cvWeighted: "
                                            "unsupported input array type");
            const unsigned* shape = histo.binContents().shapeData();
            if (!isNumpyShapeCompatible(input, shape, histo.dim()) ||
                !isNumpyShapeCompatible(leaveOneOut, shape, histo.dim()))
                throw std::invalid_argument("In npstat::AbsBandwidthCVGND::cvWeighted: "
                                            "incompatible array shape");
            double* in = (double*)PyArray_DATA((PyArrayObject*)input);
            double* leave = (double*)PyArray_DATA((PyArrayObject*)leaveOneOut);
            const ArrayND<double>& wrap = externalMemArrayND(in, shape, histo.dim());
            const ArrayND<double>& wrapL = externalMemArrayND(leave, shape, histo.dim());
            return this->operator()(histo, effectiveSampleSize, wrap, wrapL, filterUsed);
        }
#endif // SWIG
    };
}
 
#endif // NPSTAT_ABSBANDWIDTHGCV_HH_