LOrPE1DVariableDegreeCVRunner.hh

Go to the documentation of this file.

#ifndef NPSTAT_LORPE1DVARIABLEDEGREECVRUNNER_HH_
#define NPSTAT_LORPE1DVARIABLEDEGREECVRUNNER_HH_
 
/*!
// \file LOrPE1DVariableDegreeCVRunner.hh
//
// \brief Simultaneous optimization of LOrPE bandwidth and degree
//        choice by cross-validation, assuming that the golden
//        section search can succeed for both variables 
//
// Author: I. Volobouev
//
// August 2022
*/
 
#include <cmath>
#include <climits>
 
#include "npstat/nm/LinInterpolatedTable1D.hh"
#include "npstat/nm/ConstSubscriptMap.hh"
#include "npstat/nm/DualAxis.hh"
#include "npstat/nm/SimpleFunctors.hh"
 
#include "npstat/stat/LOrPE1DFixedDegreeCVRunner.hh"
 
namespace npstat {
    namespace Private {
        template<class Lorpe>
        class LOrPE1DVariableDegreeCVRunnerHelper : public Functor1<double,double>
        {
        public:
            typedef LOrPE1DCVResult Status;
            typedef ConstSubscriptMap<double,Status> StatusMap;
 
            inline LOrPE1DVariableDegreeCVRunnerHelper(
                Lorpe& lorpe, const LinInterpolatedTable1D& relativeBandwidths,
                const double referenceDegree, const double referenceBandwidth,
                const double tol)
                : lorpe_(lorpe),
                  relativeBandwidths_(relativeBandwidths),
                  referenceDegree_(referenceDegree),
                  referenceBandwidth_(referenceBandwidth),
                  tol_(tol)
            {
                assert(referenceDegree_ >= 0.0);
                assert(referenceBandwidth_ > 0.0);
                assert(tol_ > 0.0);
 
                referenceBwFactor_ = relativeBandwidths_(referenceDegree_);
                assert(referenceBwFactor_ > 0.0);
            }
 
            inline virtual ~LOrPE1DVariableDegreeCVRunnerHelper() {}
 
            inline virtual double operator()(const double& degree) const
            {
                // Try to guess a good bandwidth value for this degree
                const double bw1 = relativeBandwidths_(degree);
                const double fr = bw1/referenceBwFactor_;
                double bwGuess = referenceBandwidth_*fr;
                assert(bwGuess > 0.0);
 
                if (!statusMap_.empty())
                {
                    StatusMap::const_iterator closest = statusMap_.upper_bound(degree);
                    if (closest != statusMap_.begin())
                        --closest;
                    const double bw0 = relativeBandwidths_(closest->first);
                    const double f0 = bw1/bw0;
                    if (std::abs(std::log(f0)) < std::abs(std::log(fr)))
                    {
                        bwGuess = closest->second.bwFactor()*f0;
                        assert(bwGuess > 0.0);
                    }
                }
 
                const LOrPE1DFixedDegreeCVRunner r(degree, bwGuess, tol_);
                const Status& status = r.crossValidate(lorpe_);
                const auto insertResult = statusMap_.insert(std::make_pair(degree,status));
                assert(insertResult.second);
                return -status.cvFunction();
            }
 
            inline const StatusMap& getStatusMap() const {return statusMap_;}
 
        private:
            Lorpe& lorpe_;
            const LinInterpolatedTable1D& relativeBandwidths_;
            double referenceDegree_;
            double referenceBandwidth_;
            double referenceBwFactor_;
            double tol_;
            mutable StatusMap statusMap_;
        };
    }
 
    class LOrPE1DVariableDegreeCVRunner
    {
    public:
        inline LOrPE1DVariableDegreeCVRunner(const LinInterpolatedTable1D& i_relativeBandwidths,
                                             const double i_maxDeg,
                                             const unsigned i_nDegsInTheGrid,
                                             const double i_initialDeg,
                                             const double i_initialBwFactorGuess,
                                             const double i_initialDegStep,
                                             const double i_bwTol)
            : relativeBandwidths_(i_relativeBandwidths),
              degAxis_(UniformAxis(i_nDegsInTheGrid, 0.0, i_maxDeg)),
              initialDeg_(i_initialDeg),
              initialBwFactorGuess_(i_initialBwFactorGuess),
              initialDegStep_(i_initialDegStep),
              bwTol_(i_bwTol)
        {
            assert(i_maxDeg > 0.0);
            assert(initialDeg_ >= 0.0);
            assert(initialDeg_ <= i_maxDeg);
            assert(initialBwFactorGuess_ > 0.0);
            assert(initialDegStep_ > 0.0);
            assert(bwTol_ > 0.0);
        }
 
        template<class Lorpe>
        inline LOrPE1DCVResult crossValidate(Lorpe& lorpe) const
        {
            typedef Private::LOrPE1DVariableDegreeCVRunnerHelper<Lorpe> Helper;
            typedef LOrPE1DCVResult Result;
 
            const LOrPE1DFixedDegreeCVRunner r(
                initialDeg_, initialBwFactorGuess_, bwTol_);
            const double referenceBandwidth = r.crossValidate(lorpe).bwFactor();
            const Helper helper(lorpe, relativeBandwidths_, initialDeg_,
                                referenceBandwidth, bwTol_);
            const double gridStep = degAxis_.intervalWidth();
            const unsigned cellStep = std::round(initialDegStep_/gridStep);
            const std::pair<unsigned,double>& where = degAxis_.getInterval(initialDeg_);
            unsigned deg0 = where.first;
            if (where.second < 0.5)
                ++deg0;
 
            unsigned imin = UINT_MAX;
            double fMinusOne, fmin, fPlusOne;
 
            const MinSearchStatus1D status = goldenSectionSearchOnAGrid(
                helper, degAxis_, deg0, cellStep,
                &imin, &fMinusOne, &fmin, &fPlusOne);
            assert(status != MIN_SEARCH_FAILED);
            assert(imin < degAxis_.nCoords());
            const double degGuess = degAxis_.coordinate(imin);
            Result resGuess = helper.getStatusMap()[degGuess];
 
            if (status == MIN_SEARCH_OK)
            {
                assert(imin);
                assert(imin + 1U < degAxis_.nCoords());
 
                // Parabolic approximation in the degree variable
                double bestDegree, extremumValue;
                const bool isMinimum = parabolicExtremum(
                    degAxis_.coordinate(imin - 1U), fMinusOne,
                    degGuess,                       fmin,
                    degAxis_.coordinate(imin + 1U), fPlusOne,
                    &bestDegree, &extremumValue);
                assert(isMinimum);
 
                const double bw0 = relativeBandwidths_(degGuess);
                const double bw1 = relativeBandwidths_(bestDegree);
                const double bwtry = resGuess.bwFactor()*bw1/bw0;
                const LOrPE1DFixedDegreeCVRunner r2(bestDegree, bwtry, bwTol_);
                Result bestGuess = r2.crossValidate(lorpe);
                Result* bestResult = bestGuess.cvFunction() > resGuess.cvFunction() ?
                                     &bestGuess : &resGuess;
                bestResult->setDegreeBoundaryFlag(false);
                return *bestResult;
            }
            else
                // Optimum is on the boundary of the degree grid
                return resGuess;
        }
 
    private:
        LinInterpolatedTable1D relativeBandwidths_;
        DualAxis degAxis_;
        double initialDeg_;
        double initialBwFactorGuess_;
        double initialDegStep_;
        double bwTol_;
    };
}
 
#endif // NPSTAT_LORPE1DVARIABLEDEGREECVRUNNER_HH_