LOrPE1DFixedDegreeCVRunner.hh

#ifndef NPSTAT_LORPE1DFIXEDDEGREECVRUNNER_HH_
#define NPSTAT_LORPE1DFIXEDDEGREECVRUNNER_HH_
 
/*!
// \file LOrPE1DFixedDegreeCVScanner.hh
//
// \brief Optimization of LOrPE bandwidth choice by cross-validation
//        using the golden section search
//
// Author: I. Volobouev
//
// August 2022
*/
 
#include <cmath>
#include <cassert>
 
#include "npstat/nm/goldenSectionSearch.hh"
#include "npstat/nm/MathUtils.hh"
 
#include "npstat/stat/LOrPE1DCVResult.hh"
 
namespace npstat {
    class LOrPE1DFixedDegreeCVRunner
    {
    public:
        inline LOrPE1DFixedDegreeCVRunner(const double i_fixedDegree,
                                          const double i_bwFactorGuess,
                                          const double i_tol)
            : fixedDegree_(i_fixedDegree),
              bwFactorGuess_(i_bwFactorGuess),
              tol_(i_tol)
        {
            assert(fixedDegree_ >= 0.0);
            assert(bwFactorGuess_ > 0.0);
            assert(tol_ > 0.0);
        }
 
        inline double fixedDegree() const {return fixedDegree_;}
        inline double bwFactorGuess() const {return bwFactorGuess_;}
        inline double tol() const {return tol_;}
 
        template<class Lorpe>
        inline LOrPE1DCVResult crossValidate(Lorpe& lorpe) const
        {
            const double oldDeg = lorpe.boundFilterDegree();
            lorpe.bindFilterDegree(fixedDegree_);
 
            const double sqr2 = 1.414213562373095;
            double bwFactor = 0.0, bestFcn, divider = 1.0;
            bool bestBwFound = false;
            auto fcn = MultiplyByConst(lorpe, -1.0);
            const unsigned maxtries = 11; // divider can be increased up to
                                          // the factor 2^((maxtries-1)/2)
            for (unsigned itry=0; itry<maxtries && !bestBwFound; ++itry)
            {
                bestBwFound = goldenSectionSearchInLogSpace(
                    fcn, bwFactorGuess_/divider, tol_, &bwFactor, &bestFcn);
                divider *= sqr2;
            }
            assert(bestBwFound);
            assert(bwFactor > 0.0);
 
            // Construct parabolic approximation to the minimum
            const double fstep = 1.0 + tol_;
            const double bwMin = bwFactor/fstep;
            const double bwMax = bwFactor*fstep;
 
            double extremumCoordinate, extremumValue;
            const bool isMinimum = parabolicExtremum(
                std::log(bwMin),    fcn(bwMin),
                std::log(bwFactor), bestFcn,
                std::log(bwMax),    fcn(bwMax),
                &extremumCoordinate, &extremumValue);
            assert(isMinimum);
 
            if (oldDeg < 0.0)
                lorpe.unbindFilterDegree();
            else
                lorpe.bindFilterDegree(oldDeg);
 
            return LOrPE1DCVResult(
                fixedDegree_, std::exp(extremumCoordinate), -extremumValue,
                true, false);
        }
 
    private:
        double fixedDegree_;
        double bwFactorGuess_;
        double tol_;
    };
}
 
#endif // NPSTAT_LORPE1DFIXEDDEGREECVRUNNER_HH_