ProfLike.cc

Go to the documentation of this file.

#include "ProfLike.h"

#include <iostream>
#include <sstream>
#include <algorithm>
#include <numeric>
#include <cfloat>
#include <boost/math/special_functions/fpclassify.hpp>

#include <utl/ErrorLogger.h>

#include <TMath.h>
#include <TF1.h>

#include <Minuit2/MnMinimize.h>
#include <Minuit2/MnHesse.h>
#include <Minuit2/FCNBase.h>
#include <Minuit2/FunctionMinimum.h>

#include <Minuit2/MnUserCovariance.h>
#include <Minuit2/MnPrint.h>
#include <Minuit2/MnMigrad.h>

using std::cout;
using std::endl;


namespace MdLDFFinderAG {

  const double ProfLike::kInitError = 0.1;

  const double ProfLike::satLikeLimit = 0.001;

  const double ProfLike::kMuPoissonApprox = 0.01;


  ProfLike::ProfLike(const unsigned int n0, const utl::TraceUI& trace)
  {
    // Convert the trace to the good old vector we used so far
    std::vector<unsigned int> k0;
    for (utl::TraceUI::ConstIterator it = trace.Begin(); it != trace.End(); ++it)
      k0.push_back(*it);

    fNumberSegments = n0;

    // Sort k0 (necessary to count the number of windows with the same number of segments on)
    std::sort(k0.begin(), k0.end());

    // find unique occurrences in k0
    std::vector<unsigned int> kuniq(k0);      // kuniq holds unique occurrences in k0
    std::vector<unsigned int>::iterator it = std::unique(kuniq.begin(), kuniq.end());
    kuniq.resize(std::distance(kuniq.begin(), it));

    // Fill fSegmentsOnMulti
    for (std::vector<unsigned int>::iterator it = kuniq.begin(); it != kuniq.end(); ++it) {
      if(*it) {  // skip k=0 windows if there are more than 1 window
        const unsigned int count = std::count(k0.begin(), k0.end(), *it);
        fSegmentsOnMulti.push_back(std::make_pair(*it,count));
      }
    }

    // if fSegmentsOnMulti is empty generate a single window with 0
    if (fSegmentsOnMulti.empty())
      fSegmentsOnMulti.push_back(std::make_pair(0, 1));

    // set the maximum number of segments on and its multiplicity
    fMaxSegments = fSegmentsOnMulti.back().first;
    fMaxMultiplicity = fSegmentsOnMulti.back().second;

    /*
    DEBUGLOG("");
    // cout << "In ProfLike::ProfLike: fSegmentsOnMulti = " << endl;
    for(SegmentsOnMulti::iterator it = fSegmentsOnMulti.begin();
        it != fSegmentsOnMulti.end(); ++it) {
      // cout << (*it).first << " " << (*it).second << endl;
    }
    */
  }


  // Calculate the global minimum
  double
  ProfLike::GlobalMin()
    const
  {
    double rvalue = 0;
    for (SegmentsOnMulti::const_iterator it1 = fSegmentsOnMulti.begin();
         it1 != fSegmentsOnMulti.end(); ++it1) {
      const double k = it1->first;
      const double m = it1->second;
      if (k != fNumberSegments && k) {
        rvalue += -m*(k*log(k / fNumberSegments) + (fNumberSegments - k)*log(1-k/fNumberSegments));
      }  // else if k=fNumberSegments or k=0 GlobalMin=0
    }

    //cout << "Global minimum = " << rvalue << endl;
    return rvalue;
  }


  // Set the minuit seed
  ROOT::Minuit2::MnUserParameters
  ProfLike::SetSeed(const double mu)
    const
  {
    std::vector<double> vseed;              // seed for the minimisation

    unsigned int ksum = 0; // sum of all k's
    for (SegmentsOnMulti::const_iterator it = fSegmentsOnMulti.begin();
         it != fSegmentsOnMulti.end(); ++it) {
      ksum += it->first * it->second;
                }

    // Fill the seed vector
    for (SegmentsOnMulti::const_iterator it = fSegmentsOnMulti.begin();
         it != fSegmentsOnMulti.end(); ++it) {
        vseed.push_back((it->first) * mu/ksum);
                }

    // Trim the last element of vseed
    vseed.pop_back();

    //for (std::vector<double>::iterator it = vseed.begin(); it != vseed.end(); ++it) {
    //    cout << "In ProfLike::SetSeed(double): seed = " << *it << endl;
    //}

    return SetMnUserParameters(vseed);
  }


  // Build MnUserParameters from the seeds vector
  ROOT::Minuit2::MnUserParameters
  ProfLike::SetMnUserParameters(const std::vector<double>& vseed)
    const
  {
    // Create Minuit user parameters
    ROOT::Minuit2::MnUserParameters upar;
    std::stringstream ss;
    unsigned int i = 1;
    for (std::vector<double>::const_iterator it = vseed.begin(); it != vseed.end(); ++it) {
      ss.str("");
      ss << "mu" << i;
      upar.Add(ss.str(), *it, kInitError);
      upar.SetLowerLimit(ss.str(), 0);
      //upar.SetLimits(ss1.str(), 0, mu);    // Set limits in the parameters
      ++i;
    }

    return upar;
  }


  double
  ProfLike::operator()(const double mu)
    const
  {
    /*
    DEBUGLOG("Starting");
    cout << "In ProfLike::operator()(double): mu = " << mu << endl;
    cout << "Number of windows " << fSegmentsOnMulti.size() << endl;
    for (SegmentsOnMulti::const_iterator it1=fSegmentsOnMulti.begin() ; it1!=fSegmentsOnMulti.end(); ++it1) {
        cout << "scints on = " << (*it1).first;
        cout << ", multiplicity = " << (*it1).second << endl;
    }
    */
    double rvalue;

    // limit for the "large mu" analytical approximation of the likelihood
    // (avoids minimization issues when mu is large)
    const double mulimit = -double(fNumberSegments) * log(satLikeLimit/double(fNumberSegments));

    // calculate the mu left for the hottest window
    double mua = 0; // sum of the mle of mu over all but the hottest window
    for (SegmentsOnMulti::const_iterator it = fSegmentsOnMulti.begin();
         it != fSegmentsOnMulti.end(); ++it) { // calculate mua

      const double k = double(it->first);
      const double m = double(it->second);

      if (it->first != fMaxSegments)
        mua += -double(fNumberSegments) * log(1 - k/fNumberSegments)*m;

    }
    const double mub = (mu - mua) / fMaxMultiplicity;

    // Try different limiting cases before attempting the profile likelihood
    if (fMaxSegments == 0) {  // if there no segments on
      rvalue = mu;
    } else if (mu < kMuPoissonApprox) {  // if mu is small use Poisson
       const double mu2 = std::max(mu, DBL_EPSILON);  // avoid mu values lower than the machine precision

       //cout << "Using a Poisson likelihood" << endl;
      const double k = GetTotalSegmentsOn();
      rvalue = (mu2 - k) - k*log(mu2 / k);  // Poisson likelihood
    } else if (fSegmentsOnMulti.size() == 1) {      // if there is only 1 time window

      const unsigned int segmentsOn = fSegmentsOnMulti[0].first;
      const unsigned int multiplicity = fSegmentsOnMulti[0].second;

      double like1 = multiplicity*Like1(mu/multiplicity, segmentsOn);

//      cout << "Using 1 window likelihood. mu = " << mu << endl;
//      cout << "Like1 = " << like1 << endl;
//      cout << "GlobalMin = " << GlobalMin() << endl;

      rvalue = like1 - GlobalMin();   // minimum wrt to the global minimum

    } else if (mub > mulimit) {    // if mu is large

      // large mu approximation (the minimisation is unstable for very large mu)
      // use likelihood of the hottest window only - assumes the other as minmized

      // cast for convenience
      const double k = fMaxSegments;
      const double n = fNumberSegments;

//      cout << "Using large mu approx. mub = " << mub << endl;

      // multiplicity already considered in mub...
      double likeMaxWindow = Like1(mub, fMaxSegments);

//      cout << "rvalue before minimum substraction = " << rvalue << endl;

      double nLogLikeMin = 0;

//      cout << "fMaxSegments = " << fMaxSegments << ", fNumberSegments = ";
//      cout << fNumberSegments << endl;

      if (fMaxSegments == fNumberSegments || fMaxSegments == 0) {
        nLogLikeMin = 0;
      } else {
        nLogLikeMin = - k*log(k/n) - (n-k)*log(1-k/n);
      }

//      cout << "nLogLikeMin = " << nLogLikeMin << endl;

      rvalue = fMaxMultiplicity*(likeMaxWindow - nLogLikeMin);
      // refer to minimum use multiplicity of the hottest window

//      cout << "rvalue after offset = " << rvalue << endl;

    } else {  // calculate the profile likelihood via Minuit minimization

//      cout << "Calculating the profile likelihood using mu = " << mu << endl;

      ROOT::Minuit2::MnUserParameters upar = SetSeed(mu); // Create the initial parameters

      LikeFCN fullLike(fNumberSegments, fSegmentsOnMulti, mu);  // Create the minimization function

      ROOT::Minuit2::MnMigrad migrad(fullLike, upar); // Create the minimizer
//      DEBUGLOG("After creating the minimizer");

      gErrorIgnoreLevel = 1001;     // set the print level of Minuit (skip info messages)

      ROOT::Minuit2::FunctionMinimum min = migrad(); // minimize
//      DEBUGLOG("After running the minimizer");

//      cout << min << endl;
//      cout << "In ProfLike::operator()(double). Valid minimum = ";
//      cout << std::boolalpha << min.IsValid() << endl;

      // check that the minimization converged
      assert(min.IsValid());
      assert(min.IsValid());
      assert(min.HasValidCovariance());
      assert(!min.HesseFailed());
      assert(!min.HasReachedCallLimit());

      rvalue = min.UserState().Fval() - GlobalMin(); // minimum wrt to the global minimum

    }

    // check the return value is finite
    assert(boost::math::isfinite(rvalue));

//    DEBUGLOG("Finishing");

    return rvalue;
  }


  double
  ProfLike::Like1(const double mu, const unsigned int k)
    const
  {
    //cout << "Starting ProfLike::Like1(double mu)" << endl;

    double nlnLike = 0;

    if (mu < 1e-6) {
      nlnLike = (!k ? 0 : 1e6); // if k!=0 set -lnL to some "ridiculous" high number
    } else if (k != fNumberSegments) {
      nlnLike = mu*(1 - double(k)/fNumberSegments) - k*log(1 - exp(-mu/fNumberSegments));
    } else {  // if k==fNumberSegments
      nlnLike = -k * log(1 - exp(-mu/fNumberSegments));  // use alternative formula to avoid divergences
    }

    //cout << "mu = " << mu << ", nlnLike = " << nlnLike << endl;

    //cout << "Ending ProfLike::Like1(double mu)" << endl;

    return nlnLike;
  }


  double
  ProfLike::GetMLE()
    const
  {
    double estimatedMuons = 0;

    for (SegmentsOnMulti::const_iterator it1 = fSegmentsOnMulti.begin();
         it1 != fSegmentsOnMulti.end(); ++it1) {

      const unsigned int k = it1->first;
      const unsigned int m = it1->second;

      assert(k != fNumberSegments); // make sure there is no saturation

      double windowMuons = (-1)*double(m) * double(fNumberSegments) * log( 1 - double(k)/double(fNumberSegments) );
      estimatedMuons += windowMuons;
    }

    return estimatedMuons;
  }


  double
  ProfLike::Inverse(const double y, const double xmin, const double xmax)
    const
  {
    TF1 f1("f1", *this, xmin, xmax, 0);
    const double x = f1.GetX(y, xmin, xmax);
    return x;
  }


  double
  ProfLike::GetErrorHigh()
    const
  {
    const double muhat = GetMLE();

    // increment in the likelihood in 0.5 from the minimum to find the limit
    const double up = 0.5;

    // find xmax such as f(xmax) > up
    double xmax = muhat + 1;
    while (this->operator()(xmax) < up)
      xmax *= 2;

    const double xup = Inverse(up, muhat, xmax);
    const double xerror = xup - muhat;

    return xerror;
  }


  double
  ProfLike::GetErrorLow()
    const
  {
    const double muhat = GetMLE();
    const double xlow = GetLowLimitUnsaturated();
    const double xerror = muhat - xlow;
    return xerror;
  }


  double
  ProfLike::GetLowLimitUnsaturated()
    const
  {
    const double muhat = GetMLE();

    // increment in the likelihood in 0.5 from the minimum to find the limit
    const double up = 0.5;

    // if the event is empty
    if (!fMaxSegments)
      return 0; // the interval length is zero

    // find xmin such as f(xmin) > up
    double xmin = muhat / 2;
    while (this->operator()(xmin) < up)
      xmin /= 2;

    const double xup = Inverse(up, xmin, muhat);

    return xup;
  }


  unsigned int
  ProfLike::GetTotalSegmentsOn()
    const
  {
    unsigned int totalSegmentsOn = 0;

    for (SegmentsOnMulti::const_iterator it1 = fSegmentsOnMulti.begin();
         it1 != fSegmentsOnMulti.end(); ++it1) {

      const unsigned int k = it1->first;
      const unsigned int m = it1->second;

      totalSegmentsOn += k*m;
    }

    return totalSegmentsOn;

  }


  double
  ProfLike::GetLowLimit()
    const
  {
    // if the event is empty
    if (!fMaxSegments)
      return 0; // the interval length is zero

    // if the module is not saturated
    if (fMaxSegments < fNumberSegments)
      return GetLowLimitUnsaturated();

    // Calculate the limit in the saturated case
    // increment in the likelihood in 0.5 from the minimum to find the limit
    const double up = 0.5;
    const double xmin = 0.1;
    const double xmaxSeed = 1;

    // find xmax such as f(xmin) > up
    double xmax = xmaxSeed;
    while (this->operator()(xmax) > up)
      xmax *= 2;

    const double limit = Inverse(up, xmin, xmax);

    return limit;
  }

}