solveForLOrPEWeights.hh

Go to the documentation of this file.

#ifndef NPSTAT_SOLVEFORLORPEWEIGHTS_HH_
#define NPSTAT_SOLVEFORLORPEWEIGHTS_HH_
 
/*!
// \file solveForLOrPEWeights.hh
//
// \brief Solver for the non-linear Fredholm equation in the semiparametric
//        mixture model with unbinned LOrPE
//
// Author: I. Volobouev
//
// July 2022
*/
 
#include <iosfwd>
#include <vector>
#include <cassert>
#include <utility>
#include <algorithm>
 
#include "npstat/stat/LOrPE1DCVResult.hh"
#include "npstat/stat/AbsDistribution1D.hh"
 
namespace npstat {
    class LOrPEWeightsResult
    {
    public:
        inline LOrPEWeightsResult(const double i_alpha,
                                  const bool i_converged,
                                  const unsigned i_nIterations,
                                  const LOrPE1DCVResult& cvResult)
            : cvResult_(cvResult),
              alpha_(i_alpha),
              nIterations_(i_nIterations),
              converged_(i_converged)
        {
            if (converged_)
                assert(cvResult_.isValid());
        }
 
        inline double filterDegree() const {return cvResult_.filterDegree();}
        inline double bwFactor() const {return cvResult_.bwFactor();}
        inline double cvFunction() const {return cvResult_.cvFunction();}
        inline bool isOnDegreeBoundary() const {return cvResult_.isOnDegreeBoundary();}
        inline bool isOnBandwidthBoundary() const {return cvResult_.isOnBandwidthBoundary();}
 
        inline double alpha() const {return alpha_;}
        inline unsigned nIterations() const {return nIterations_;}
        inline bool converged() const {return converged_;}
 
        void print(std::ostream& os) const;
 
    private:
        LOrPE1DCVResult cvResult_;
        double alpha_;
        unsigned nIterations_;
        bool converged_;
    };
 
 
    template<class Lorpe, class CvRunner, class ConvergenceCalculator>
    inline LOrPEWeightsResult solveForLOrPEWeights(
        Lorpe& lorpe, const CvRunner& cvRunner,
        const AbsDistribution1D& signal, const double alpha,
        const ConvergenceCalculator& conv, const unsigned maxIterations)
    {
        const auto& sample = lorpe.coords();
        const unsigned long sampleSize = sample.size();
 
        std::vector<double> signalDensity(sampleSize);
        for (unsigned long i=0; i<sampleSize; ++i)
            signalDensity[i] = signal.density(sample[i].first);
 
        std::vector<double> weightBuffer(sampleSize*2);
        double* previousWeights = &weightBuffer[0];
        double* currentWeights = previousWeights + sampleSize;
        for (unsigned long i=0; i<sampleSize; ++i)
            previousWeights[i] = sample[i].second;
 
        std::vector<double> bgDensityBuffer(sampleSize*2);
        double* previousBg = &bgDensityBuffer[0];
        double* currentBg = previousBg + sampleSize;
 
        // Run the initial cross-validation
        LOrPE1DCVResult cvResult = cvRunner.crossValidate(lorpe);
 
        {
            // Create the initial background estimate
            // at the sample points
            const double filterDegree = cvResult.filterDegree();
            const double bwFactor = cvResult.bwFactor();
 
            for (unsigned long i=0; i<sampleSize; ++i)
            {
                previousBg[i] = lorpe.density(
                    sample[i].first, filterDegree, bwFactor);
                assert(previousBg[i] >= 0.0);
            }
        }
 
        bool converged = false;
        unsigned iter = 0;
        for (; iter<maxIterations && !converged; ++iter)
        {
            // Update the weights
            for (unsigned long i=0; i<sampleSize; ++i)
            {
                if (previousBg[i] > 0.0)
                {
                    const double d = alpha*signalDensity[i] +
                                     (1.0 - alpha)*previousBg[i];
                    assert(d > 0.0);
                    currentWeights[i] = previousBg[i]/d;
                }
                else
                    currentWeights[i] = 0.0;
            }
            lorpe.setWeights(currentWeights, sampleSize);
 
            // Run cross-validation with the new weights
            cvResult = cvRunner.crossValidate(lorpe);
            const double filterDegree = cvResult.filterDegree();
            const double bwFactor = cvResult.bwFactor();
 
            // Calculate background density
            for (unsigned long i=0; i<sampleSize; ++i)
            {
                currentBg[i] = lorpe.density(
                    sample[i].first, filterDegree, bwFactor);
                assert(currentBg[i] >= 0.0);
            }
 
            // Check for convergence
            converged = conv.converged(lorpe,
                                       previousWeights, currentWeights,
                                       previousBg, currentBg);
 
            // Rotate buffers
            std::swap(previousWeights, currentWeights);
            std::swap(previousBg, currentBg);
        }
 
        return LOrPEWeightsResult(alpha, converged, iter, cvResult);
    }
}
 
std::ostream& operator<<(std::ostream& os, const npstat::LOrPEWeightsResult& res);
 
#endif // NPSTAT_SOLVEFORLORPEWEIGHTS_HH_