AbsUnfold1D.hh

Go to the documentation of this file.

#ifndef NPSTAT_ABSUNFOLD1D_HH_
#define NPSTAT_ABSUNFOLD1D_HH_
 
/*!
// \file AbsUnfold1D.hh
//
// \brief Interface definition for 1-d unfolding algorithms
//
// Author: I. Volobouev
//
// February 2014
*/
 
#include <vector>
#include <utility>
#include <cassert>
#include <stdexcept>
 
#include "npstat/nm/Matrix.hh"
 
#ifdef SWIG
#include "npstat/wrap/NumpyTypecode.hh"
#endif // SWIG
 
namespace npstat {
    class LocalPolyFilter1D;
 
    class AbsUnfold1D
    {
    public:
        /**
        // Number of rows of the response matrix should be equal
        // to the number of discretization intervals in the 
        // space of observations. Number of columns should be equal
        // to the number of intervals in the original (unfolded) space.
        */
        explicit AbsUnfold1D(const Matrix<double>& responseMatrix);
 
        inline virtual ~AbsUnfold1D() {}
 
        inline const Matrix<double>& responseMatrix() const
            {return responseMatrix_;}
 
        inline const Matrix<double>& transposedResponseMatrix() const
            {return responseMatrixT_;}
 
        inline const std::vector<double>& efficiency() const
            {return efficiency_;}
 
        //@{
        /**
        // The first element of the pair is the entropic nDoF and
        // the second is the nDoF based on the trace
        */
        std::pair<double, double> smoothingNDoF() const;
        std::pair<double, double> responseNDoF() const;
        std::pair<double, double> smoothedResponseNDoF() const;
        //@}
 
        /**
        // Effective number of degrees of freedom for the smearing model
        // (for use in AIC, etc). Will call "unfold" internally, so can't
        // be a const method.
        */
        std::pair<double, double> modelNDoF();
 
        /** Set the initial approximation to the unfolded solution */
        virtual void setInitialApproximation(
            const double* approx, unsigned lenApprox);
 
        /** Clear the initial approximation to the unfolded solution */
        virtual void clearInitialApproximation();
 
        /** Return the initial approximation to the unfolded solution */
        virtual const std::vector<double>& getInitialApproximation() const;
 
        /**
        // Set the smoothing filter used. The filter will not be copied.
        // It is the responsibility of the user of this class to ensure
        // that the lifetime of the filter exceeds the lifetime of
        // this object.
        */
        virtual void setFilter(const LocalPolyFilter1D* f);
 
        /** Retrieve the smoothing filter used */
        virtual const LocalPolyFilter1D* getFilter(bool throwIfNull=false) const;
 
        /** Switch between using filtering or convolution */
        inline virtual void useConvolutions(const bool b) {useConvolutions_ = b;}
 
        /** Check if the filter should use "filter" or "convolve" method */
        inline bool usingConvolutions() const {return useConvolutions_;}
 
        /**
        // Method to be implemented by derived classes. "lenObserved"
        // and "lenUnfolded" should be consistent with the dimensionality
        // of the response matrix provided in the constructor. The
        // "observationCovarianceMatrix" pointer to the covariance matrix
        // of observations can be NULL in which case this covariance
        // matrix will be evaluated internally (typically, assuming
        // either Poisson or multinomial distributions for the observed
        // counts). "unfoldedCovarianceMatrix" pointer, if not NULL,
        // should be used to return the covariance matrix of unfolded
        // values with dimensions (lenUnfolded x lenUnfolded). This
        // function should return "true" on success and "false" on failure.
        */
        virtual bool unfold(const double* observed, unsigned lenObserved,
                            const Matrix<double>* observationCovarianceMatrix,
                            double* unfolded, unsigned lenUnfolded,
                            Matrix<double>* unfoldedCovarianceMatrix) = 0;
 
        /** 
        // An L1-like relative distance between two unnormalized
        // density-like arrays. Can be used as a convergence measure.
        */
        static double probDelta(const double* prev, const double* next,
                                const unsigned len);
 
        /** 
        // Procedure for building covariance matrix of observations
        // according to either Poisson or multinomial distribution
        */
        static Matrix<double> observationCovariance(
            const double* yhat, unsigned lenObserved, bool isMultinomial);
 
    protected:
        // Build uniform initial approximation to the unfolded solution.
        // This is useful if the initial approximation was not set explicitly.
        void buildUniformInitialApproximation(const double* observed,
                                              unsigned lenObserved,
                                              std::vector<double>* result) const;
 
        // The following function will throw the std::invalid_argument
        // exception if the dimensions are incompatible with those of the
        // response matrix
        void validateDimensions(unsigned lenObserved, unsigned lenUnfold) const;
 
        // The following function will throw the std::invalid_argument
        // exception if the input array has negative entries or if all
        // values are 0.
        static void validateDensity(const double* observ,
                                    unsigned lenObserved);
    private:
        // Disable copy and assignment because using the same
        // filter by two different objects of this type is not
        // necessarily cool
        AbsUnfold1D();
        AbsUnfold1D(const AbsUnfold1D&);
        AbsUnfold1D& operator=(const AbsUnfold1D&);
 
        const Matrix<double> responseMatrix_;
        const Matrix<double> responseMatrixT_;
        std::vector<double> efficiency_;
        std::vector<double> initialApproximation_;
        const LocalPolyFilter1D* filt_;
        bool useConvolutions_;
 
#ifdef SWIG
    public:
        inline const LocalPolyFilter1D* getFilter2() const
        {
            return this->getFilter(true);
        }
        static inline double probDelta2(const double* prev, unsigned prevlen,
                                        const double* next, unsigned nextlen)
        {
            if (prevlen != nextlen) throw std::runtime_error(
                "In npstat::AbsUnfold1D::probDelta2: incompatible inputs");
            return probDelta(prev, next, nextlen);
        }
        inline PyObject* unfold2(const double* observed, unsigned lenObserved,
                                 Matrix<double>* covMatOut)
        {
            assert(covMatOut);
            const int typenum = NumpyTypecode<double>::code;
            npy_intp sh = responseMatrix_.nColumns();
            PyObject* array = PyArray_SimpleNew(1, &sh, typenum);
            double* result = (double*)PyArray_DATA((PyArrayObject*)array);
            try
            {
                const bool status = this->unfold(observed, lenObserved, 0,
                                                 result, sh, covMatOut);
                if (!status) throw std::runtime_error(
                    "In npstat::AbsUnfold1D::unfold2: "
                    "expectation-maximization iterations failed to converge");
            }
            catch (const std::exception& e)
            {
                Py_DECREF(array);
                throw;
            }
            return array;
        }
        inline PyObject* unfold2(const double* observed, unsigned lenObserved,
                                 const Matrix<double>& observationCovMat,
                                 Matrix<double>* covMatOut)
        {
            assert(covMatOut);
            const int typenum = NumpyTypecode<double>::code;
            npy_intp sh = responseMatrix_.nColumns();
            PyObject* array = PyArray_SimpleNew(1, &sh, typenum);
            double* result = (double*)PyArray_DATA((PyArrayObject*)array);
            try
            {
                const bool status = this->unfold(observed, lenObserved,
                                                 &observationCovMat,
                                                 result, sh, covMatOut);
                if (!status) throw std::runtime_error(
                    "In npstat::AbsUnfold1D::unfold2: "
                    "expectation-maximization iterations failed to converge");
            }
            catch (const std::exception& e)
            {
                Py_DECREF(array);
                throw;
            }
            return array;
        }
        inline PyObject* unfoldNoCov(const double* observed, unsigned lenObserved)
        {
            const int typenum = NumpyTypecode<double>::code;
            npy_intp sh = responseMatrix_.nColumns();
            PyObject* array = PyArray_SimpleNew(1, &sh, typenum);
            double* result = (double*)PyArray_DATA((PyArrayObject*)array);
            try
            {
                const bool status = this->unfold(observed, lenObserved, 0,
                                                 result, sh, 0);
                if (!status) throw std::runtime_error(
                    "In npstat::AbsUnfold1D::unfoldNoCov: "
                    "expectation-maximization iterations failed to converge");
            }
            catch (const std::exception& e)
            {
                Py_DECREF(array);
                throw;
            }
            return array;
        }
        inline PyObject* unfoldNoCov(const double* observed, unsigned lenObserved,
                                     const Matrix<double>& observationCovMat)
        {
            const int typenum = NumpyTypecode<double>::code;
            npy_intp sh = responseMatrix_.nColumns();
            PyObject* array = PyArray_SimpleNew(1, &sh, typenum);
            double* result = (double*)PyArray_DATA((PyArrayObject*)array);
            try
            {
                const bool status = this->unfold(observed, lenObserved,
                                                 &observationCovMat,
                                                 result, sh, 0);
                if (!status) throw std::runtime_error(
                    "In npstat::AbsUnfold1D::unfoldNoCov: "
                    "expectation-maximization iterations failed to converge");
            }
            catch (const std::exception& e)
            {
                Py_DECREF(array);
                throw;
            }
            return array;
        }
        inline PyObject* fold(const double* unfolded, unsigned lenUnfolded) const
        {
            if (lenUnfolded != responseMatrix_.nColumns()) throw std::invalid_argument(
                "In npstat::AbsUnfold1D::fold: incompatible input array size");
            const int typenum = NumpyTypecode<double>::code;
            npy_intp sh = responseMatrix_.nRows();
            PyObject* array = PyArray_SimpleNew(1, &sh, typenum);
            double* result = (double*)PyArray_DATA((PyArrayObject*)array);
            responseMatrix_.timesVector(unfolded, lenUnfolded, result, sh);
            return array;
        }
#endif // SWIG
    };
}
 
#endif // NPSTAT_ABSUNFOLD1D_HH_