LOrPE1DVariableDegreeCVPicker.hh

Go to the documentation of this file.

#ifndef NPSTAT_LORPE1DVARIABLEDEGREECVPICKER_HH_
#define NPSTAT_LORPE1DVARIABLEDEGREECVPICKER_HH_
 
/*!
// \file LOrPE1DVariableDegreeCVPicker.hh
//
// \brief Simultaneous optimization of LOrPE bandwidth and degree
//        choice by cross-validation
//
// Author: I. Volobouev
//
// August 2022
*/
 
#include "npstat/nm/LinearMapper1d.hh"
 
#include "npstat/stat/LOrPE1DFixedDegreeCVPicker.hh"
 
namespace npstat {
    class LOrPE1DVariableDegreeCVPicker;
 
    namespace Private {
        template<class Lorpe>
        class LOrPE1DVariableDegreeCVPickerHelper : public Functor1<double,double>
        {
        public:
            typedef std::map<double,LOrPE1DCVResult> StatusMap;
 
            inline LOrPE1DVariableDegreeCVPickerHelper(
                Lorpe& lorpe, const LOrPE1DVariableDegreeCVPicker& picker)
                : lorpe_(lorpe),
                  picker_(picker)
            {}
 
            inline virtual ~LOrPE1DVariableDegreeCVPickerHelper() {}
 
            virtual double operator()(const double& degree) const;
 
            // This method returns an invalid result if the degree key is not memoized
            inline LOrPE1DCVResult getMemoized(const double degree) const
            {
                const auto it = statusMap_.find(degree);
                if (it == statusMap_.end())
                    return LOrPE1DCVResult();
                else
                    return it->second;
            }
 
        private:
            Lorpe& lorpe_;
            const LOrPE1DVariableDegreeCVPicker& picker_;
            mutable StatusMap statusMap_;
        };
    }
 
    class LOrPE1DVariableDegreeCVPicker
    {
    public:
        /**
        // LOrPE1DVariableDegreeCVPicker constructor arguments are:
        //
        // i_maxDeg                     -- Maximum filter degree to consider
        //                                 in the optimization procedure.
        //
        // minBwFactors                 -- Minimum bandwidth factors to
        //                                 consider for each filter degree.
        //                                 The number of points in the degree
        //                                 grid will be set to the size
        //                                 of this vector. For example,
        //                                 if i_maxDeg = 2 and the size
        //                                 of minBwFactors is 5, the following
        //                                 five degrees will be considered:
        //                                 0.0, 0.5, 1.0, 1.5, 2.0. The
        //                                 first element of "minBwFactors"
        //                                 will be used for degree 0, the
        //                                 second for degree 0.5, etc. The
        //                                 size of "minBwFactors" vector
        //                                 must be at least 2.
        //
        // maxBwFactors                 -- Maximum bandwidth factors to
        //                                 consider for each degree. The size
        //                                 of this vector should be equal to
        //                                 the size of "minBwFactors".
        //
        // startingBwFactors            -- Bandwidth factors to use to start
        //                                 searches (for each filter degree).
        //                                 The size of this vector should be
        //                                 equal to the size of "minBwFactors".
        //
        // nBwFactors                   -- The size of the grid of bandwidth
        //                                 factors. For each degree, this grid
        //                                 will be generated by
        //                                 "EquidistantInLogSpace" class with
        //                                 the minimum and maximum taken from
        //                                 the corresponding elements of
        //                                 "minBwFactors" and "maxBwFactors"
        //                                 vectors. This parameter must be
        //                                 at least 2.
        //
        // initialStepSizeInFactorCells -- Initial step size in the units of
        //                                 bandwidth cell width, for the
        //                                 optimization over bandwidth factors.
        //
        // initialDeg                   -- Initial guess for the best filter
        //                                 degree.
        //
        // initialStepSizeInDegreeCells -- Initial step size in the units of
        //                                 filter degree cell width, for the
        //                                 optimization over filter degrees.
        //                                 The width of the filter degree cell
        //                                 is just
        //                                 i_maxDeg/(minBwFactors.size() - 1).
        */
        inline LOrPE1DVariableDegreeCVPicker(const double i_maxDeg,
                                             const std::vector<double>& minBwFactors,
                                             const std::vector<double>& maxBwFactors,
                                             const std::vector<double>& startingBwFactors,
                                             const unsigned nBwFactors,
                                             const unsigned initialStepSizeInFactorCells,
                                             const double initialDeg,
                                             const unsigned initialStepSizeInDegreeCells)
            : degAxis_(UniformAxis(minBwFactors.size(), 0.0, i_maxDeg)),
              minBwFactorForEachDeg_(minBwFactors),
              maxBwFactorForEachDeg_(maxBwFactors),
              startingBwFactorForEachDeg_(startingBwFactors),
              nBwFactors_(nBwFactors),
              initialStepSizeInFactorsGridCells_(initialStepSizeInFactorCells),
              initialStepSizeInDegreeGridCells_(initialStepSizeInDegreeCells)
        {
            const unsigned nDegsInTheGrid = minBwFactorForEachDeg_.size();
 
            assert(i_maxDeg > 0.0);
            assert(nDegsInTheGrid > 1U);
            assert(maxBwFactorForEachDeg_.size() == nDegsInTheGrid);
            assert(startingBwFactorForEachDeg_.size() == nDegsInTheGrid);
 
            for (unsigned i=0; i<nDegsInTheGrid; ++i)
            {
                assert(minBwFactorForEachDeg_[i] > 0.0);
                assert(maxBwFactorForEachDeg_[i] > minBwFactorForEachDeg_[i]);
                assert(startingBwFactorForEachDeg_[i] > 0.0);
            }
            assert(nBwFactors_ > 1U);
 
            if (initialDeg <= 0.0)
                initialDegCell_ = 0U;
            else if (initialDeg >= i_maxDeg)
                initialDegCell_ = nDegsInTheGrid - 1U;
            else
            {
                const std::pair<unsigned,double>& cell = degAxis_.getInterval(initialDeg);
                initialDegCell_ = cell.first;
                if (cell.second < 0.5)
                    ++initialDegCell_;
            }
 
            if (!initialStepSizeInDegreeGridCells_)
                initialStepSizeInDegreeGridCells_ = 1U;
            if (!initialStepSizeInFactorsGridCells_)
                initialStepSizeInFactorsGridCells_ = 1U;
        }
 
        template<class Lorpe>
        inline LOrPE1DCVResult crossValidate(Lorpe& lorpe) const
        {
            typedef Private::LOrPE1DVariableDegreeCVPickerHelper<Lorpe> Helper;
            typedef LOrPE1DCVResult Result;
 
            const unsigned nscan = degAxis_.nCoords();
            const Helper helper(lorpe, *this);
            unsigned imin = UINT_MAX;
            double fMinusOne, fmin, fPlusOne;
 
            const MinSearchStatus1D status = goldenSectionSearchOnAGrid(
                helper, degAxis_, initialDegCell_, initialStepSizeInDegreeGridCells_,
                &imin, &fMinusOne, &fmin, &fPlusOne);
 
            if (status == MIN_SEARCH_FAILED)
            {
                std::vector<double> buffer(2*nscan);
                double* degrees = &buffer[0];
                double* cvValues = degrees + nscan;
 
                for (unsigned i=0; i<nscan; ++i)
                {
                    const double deg = degAxis_.coordinate(i);
                    degrees[i] = deg;
                    const LOrPE1DCVResult& searched = helper.getMemoized(deg);
                    if (searched.isValid())
                        cvValues[i] = searched.cvFunction();
                    else
                        cvValues[i] = -helper(deg);
                }
 
                const ScanExtremum1D& scanMax = findScanMaximum1D(degrees, nscan,
                                                                  cvValues, nscan);
                LOrPE1DCVResult searched = helper.getMemoized(scanMax.location());
                if (scanMax.isOnTheBoundary())
                    assert(searched.isValid());
                else if (searched.isValid())
                    searched.setDegreeBoundaryFlag(false);
                else
                    // Run actuall cross-validation at the guessed degree
                    searched = cvArbitraryDegree(lorpe, scanMax.location());
                return searched;
            }
            else if (status == MIN_SEARCH_OK)
            {
                assert(imin);
                assert(imin + 1U < nscan);
 
                const double degGuess = degAxis_.coordinate(imin);
                Result resGuess = helper.getMemoized(degGuess);
                assert(resGuess.isValid());
                resGuess.setDegreeBoundaryFlag(false);
 
                // 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);
 
                // Run actuall cross-validation at the guessed degree
                const Result bestGuess = cvArbitraryDegree(lorpe, bestDegree);
                return bestGuess.cvFunction() > resGuess.cvFunction() ?
                       bestGuess : resGuess;
            }
            else
            {
                // Optimum is on the boundary of the degree grid
                assert(imin < nscan);
                const double deg = degAxis_.coordinate(imin);
                const LOrPE1DCVResult& searched = helper.getMemoized(deg);
                assert(searched.isValid());
                return searched;
            }
        }
 
    private:
        template<class Lorpe>
        friend class Private::LOrPE1DVariableDegreeCVPickerHelper;
 
        template<class Lorpe>
        inline LOrPE1DCVResult cvArbitraryDegree(Lorpe& lorpe, const double bestDegree) const
        {
            const unsigned nscan = degAxis_.nCoords();
 
            double minBwFactor, maxBwFactor, firstBwFactorToTry;
            const unsigned ibest = degAxis_.getInterval(bestDegree).first;
            const unsigned inext = ibest + 1U;
            if (inext < nscan)
            {
                const double xbest = degAxis_.coordinate(ibest);
                const double xnext = degAxis_.coordinate(inext);
 
                // Perform logarithmic interpolation of the min/max/firstToTry
                // bandwidth factors to the given degree which no longer
                // has to be on the degree grid
                const LinearMapper1d m0(xbest, std::log(minBwFactorForEachDeg_[ibest]),
                                        xnext, std::log(minBwFactorForEachDeg_[inext]));
                minBwFactor = std::exp(m0(bestDegree));
 
                const LinearMapper1d m1(xbest, std::log(maxBwFactorForEachDeg_[ibest]),
                                        xnext, std::log(maxBwFactorForEachDeg_[inext]));
                maxBwFactor = std::exp(m1(bestDegree));
 
                const LinearMapper1d m2(xbest, std::log(startingBwFactorForEachDeg_[ibest]),
                                        xnext, std::log(startingBwFactorForEachDeg_[inext]));
                firstBwFactorToTry = std::exp(m2(bestDegree));
            }
            else
            {
                minBwFactor = minBwFactorForEachDeg_[ibest];
                maxBwFactor = maxBwFactorForEachDeg_[ibest];
                firstBwFactorToTry = startingBwFactorForEachDeg_[ibest];
            }
 
            const LOrPE1DFixedDegreeCVPicker r2(bestDegree,
                    minBwFactor, maxBwFactor, nBwFactors_,
                    firstBwFactorToTry, initialStepSizeInFactorsGridCells_);
            LOrPE1DCVResult cv = r2.crossValidate(lorpe);
            assert(cv.isValid());
            cv.setDegreeBoundaryFlag(false);
            return cv;
        }
 
        DualAxis degAxis_;
        std::vector<double> minBwFactorForEachDeg_;
        std::vector<double> maxBwFactorForEachDeg_;
        std::vector<double> startingBwFactorForEachDeg_;
        unsigned nBwFactors_;
        unsigned initialStepSizeInFactorsGridCells_;
        unsigned initialDegCell_;
        unsigned initialStepSizeInDegreeGridCells_;
    };
 
    namespace Private {
        template<class Lorpe>
        inline double
        LOrPE1DVariableDegreeCVPickerHelper<Lorpe>::operator()(const double& degree) const
        {
            // Check that the degree argument is actually on
            // the degree grid. Figure out the bandwidth range
            // and the first bandwidth to try for this degree.
            const double eps = 1.e-5;
            const std::pair<unsigned,double>& cell = picker_.degAxis_.getInterval(degree);
            assert(cell.second < eps || cell.second > 1.0 - eps);
            unsigned ideg = cell.first;
            if (cell.second < eps)
                ++ideg;
            assert(ideg < picker_.degAxis_.nCoords());
 
            const double minBwFactor = picker_.minBwFactorForEachDeg_[ideg];
            const double maxBwFactor = picker_.maxBwFactorForEachDeg_[ideg];
            const double firstBwFactorToTry = picker_.startingBwFactorForEachDeg_[ideg];
 
            // Now, find the optimal bandwidth for this degree
            const LOrPE1DFixedDegreeCVPicker r(
                degree, minBwFactor, maxBwFactor, picker_.nBwFactors_,
                firstBwFactorToTry, picker_.initialStepSizeInFactorsGridCells_);
            const LOrPE1DCVResult& status = r.crossValidate(lorpe_);
            assert(status.isValid());
 
            // Remember the optimal result for this degree
            const auto insertResult = statusMap_.insert(std::make_pair(degree,status));
            assert(insertResult.second);
            return -status.cvFunction();
        }
    }
}
 
#endif // NPSTAT_LORPE1DVARIABLEDEGREECVPICKER_HH_