LOrPE1DFixedDegreeCVPicker.hh

Go to the documentation of this file.

#ifndef NPSTAT_LORPE1DFIXEDDEGREECVPICKER_HH_
#define NPSTAT_LORPE1DFIXEDDEGREECVPICKER_HH_
 
/*!
// \file LOrPE1DFixedDegreeCVPicker.hh
//
// \brief Optimization of LOrPE bandwidth choice by cross-validation
//        with a hybrid algorithm combining golden section search and
//        a grid scan
//
// Author: I. Volobouev
//
// August 2022
*/
 
#include <cmath>
#include <cassert>
#include <climits>
#include <stdexcept>
#include <map>
 
#include "npstat/nm/goldenSectionSearch.hh"
#include "npstat/nm/MathUtils.hh"
#include "npstat/nm/ScanExtremum1D.hh"
#include "npstat/nm/EquidistantSequence.hh"
#include "npstat/nm/DualAxis.hh"
#include "npstat/nm/SimpleFunctors.hh"
 
#include "npstat/stat/LOrPE1DCVResult.hh"
 
namespace npstat {
    namespace Private {
        template<class Lorpe>
        class LOrPE1DFixedDegreeCVPickerHelper : public Functor1<double,double>
        {
        public:
            inline LOrPE1DFixedDegreeCVPickerHelper(Lorpe& lorpe, const double fixedDegree)
                : lorpe_(lorpe),
                  fixedDegree_(fixedDegree)
            {
                assert(fixedDegree_ >= 0.0);
            }
 
            inline virtual ~LOrPE1DFixedDegreeCVPickerHelper() {}
 
            inline virtual double operator()(const double& bwFactor) const
            {
                const double cv = lorpe_(fixedDegree_, bwFactor);
                const auto insertResult = cvValues_.insert(std::make_pair(bwFactor,cv));
                assert(insertResult.second);
                return -cv;
            }
 
            inline std::pair<double,bool> getMemoized(const double bwFactor) const
            {
                const auto it = cvValues_.find(bwFactor);
                if (it == cvValues_.end())
                    return std::make_pair(0.0, false);
                else
                    return std::make_pair(it->second, true);
            }
 
        private:
            Lorpe& lorpe_;
            double fixedDegree_;
            mutable std::map<double,double> cvValues_;
        };
    }
 
    class LOrPE1DFixedDegreeCVPicker
    {
    public:
        inline LOrPE1DFixedDegreeCVPicker(const double degree,
                                          const double minBwFactor,
                                          const double maxBwFactor,
                                          const unsigned nBwFactors,
                                          const double firstBwFactorToTry,
                                          const unsigned initialStepSizeInFactorsGridCells)
            : bandwidthFactors_(EquidistantInLogSpace(minBwFactor, maxBwFactor, nBwFactors)),
              fixedDegree_(degree),
              initialStep_(initialStepSizeInFactorsGridCells)
        {
            if (fixedDegree_ < 0.0) throw std::invalid_argument(
                "In npstat::LOrPE1DFixedDegreeCVPicker constructor: "
                "degree argument must be non-negative");
            if (minBwFactor >= maxBwFactor) throw std::invalid_argument(
                "In npstat::LOrPE1DFixedDegreeCVPicker constructor: "
                "minimum factor must be less than maximum");
            i0_ = bandwidthFactors_.getInterval(firstBwFactorToTry).first;
            assert(i0_ < nBwFactors);
            if (!initialStep_)
                initialStep_ = 1U;
        }
 
        inline unsigned nBwFactors() const {return bandwidthFactors_.nCoords();}
        inline double getBwFactor(const unsigned i) const {return bandwidthFactors_.coordinate(i);}
        inline double searchStart() const {return bandwidthFactors_.coordinate(i0_);}
 
        template<class Lorpe>
        inline LOrPE1DCVResult crossValidate(Lorpe& lorpe) const
        {
            typedef Private::LOrPE1DFixedDegreeCVPickerHelper<Lorpe> Helper;
 
            const unsigned nscan = bandwidthFactors_.nCoords();
            const Helper helper(lorpe, fixedDegree_);
            unsigned imin = UINT_MAX;
            double fMinusOne, fmin, fPlusOne;
 
            const MinSearchStatus1D status = goldenSectionSearchOnAGrid(
                helper, bandwidthFactors_, i0_, initialStep_,
                &imin, &fMinusOne, &fmin, &fPlusOne);
 
            if (status == MIN_SEARCH_FAILED)
            {
                // This most likely means we have found a local maximum
                // instead of a minimum. Switch to a simple scan instead
                // to find the global minimum.
                std::vector<double> buffer(2*nscan);
                double* bwLogs = &buffer[0];
                double* cvValues = bwLogs + nscan;
 
                for (unsigned i=0; i<nscan; ++i)
                {
                    const double bwFactor = bandwidthFactors_.coordinate(i);
                    bwLogs[i] = std::log(bwFactor);
                    const std::pair<double,bool>& searched = helper.getMemoized(bwFactor);
                    if (searched.second)
                        cvValues[i] = searched.first;
                    else
                        cvValues[i] = lorpe(fixedDegree_, bwFactor);
                }
 
                const ScanExtremum1D& scanMax = findScanMaximum1D(bwLogs, nscan,
                                                                  cvValues, nscan);
                return LOrPE1DCVResult(
                    fixedDegree_, std::exp(scanMax.location()), scanMax.value(),
                    true, scanMax.isOnTheBoundary());
            }
            else if (status == MIN_SEARCH_OK)
            {
                assert(imin);
                assert(imin + 1U < nscan);
 
                const double bwMin = bandwidthFactors_.coordinate(imin-1U);
                const double bwFactor = bandwidthFactors_.coordinate(imin);
                const double bwMax =  bandwidthFactors_.coordinate(imin+1U);
 
                // Build parabolic approximation in the log space
                double extremumCoordinate, extremumValue;
                const bool isMinimum = parabolicExtremum(
                    std::log(bwMin),    fMinusOne,
                    std::log(bwFactor), fmin,
                    std::log(bwMax),    fPlusOne,
                    &extremumCoordinate, &extremumValue);
                assert(isMinimum);
 
                return LOrPE1DCVResult(
                    fixedDegree_, std::exp(extremumCoordinate), -extremumValue,
                    true, false);
            }
            else
            {
                // The minimum is on the grid boundary
                assert(imin < nscan);
                const double bwFactor = bandwidthFactors_.coordinate(imin);
                const std::pair<double,bool>& searched = helper.getMemoized(bwFactor);
                assert(searched.second);
                return LOrPE1DCVResult(
                    fixedDegree_, bwFactor, searched.first, true, true);
            }
        }
 
    private:
        DualAxis bandwidthFactors_;
        double fixedDegree_;
        unsigned i0_;
        unsigned initialStep_;
    };
}
 
#endif // NPSTAT_LORPE1DFIXEDDEGREECVPICKER_HH_