ProcLikelihood.cc

Go to the documentation of this file.

#include <algorithm>
#include <iostream>
#include <numeric>
#include <iomanip>
#include <cstring>
#include <vector>
#include <string>
#include <memory>
#include <map>

#include <xercesc/dom/DOM.hpp>

#include <TH1.h>

#include "FWCore/Utilities/interface/Exception.h"

#include "PhysicsTools/MVAComputer/interface/AtomicId.h"

#include "PhysicsTools/MVATrainer/interface/XMLSimpleStr.h"
#include "PhysicsTools/MVATrainer/interface/XMLUniStr.h"
#include "PhysicsTools/MVATrainer/interface/XMLDocument.h"
#include "PhysicsTools/MVATrainer/interface/MVATrainer.h"
#include "PhysicsTools/MVATrainer/interface/TrainProcessor.h"

XERCES_CPP_NAMESPACE_USE

using namespace PhysicsTools;

namespace { // anonymous

class ProcLikelihood : public TrainProcessor {
    public:
    typedef TrainProcessor::Registry<ProcLikelihood>::Type Registry;

    ProcLikelihood(const char *name, const AtomicId *id,
                   MVATrainer *trainer);
    ~ProcLikelihood() override;

    void configure(DOMElement *elem) override;
    Calibration::VarProcessor *getCalibration() const override;

    void trainBegin() override;
    void trainData(const std::vector<double> *values,
                           bool target, double weight) override;
    void trainEnd() override;

    bool load() override;
    void save() override;

    struct PDF {
        std::vector<double>     distr;
        Calibration::HistogramD::Range  range;
    };

    private:
    enum Iteration {
        ITER_EMPTY,
        ITER_RANGE,
        ITER_FILL,
        ITER_DONE
    };

    struct SigBkg {
        PDF     signal;
        PDF     background;
        unsigned int    smooth;
        Iteration   iteration;
    };

    std::vector<SigBkg> pdfs;
    std::vector<double> sigSum;
    std::vector<double> bkgSum;
    std::vector<double> bias;
    int         categoryIdx;
    bool            logOutput;
    bool            individual;
    bool            neverUndefined;
    bool            keepEmpty;
    unsigned int        nCategories;
    bool            doCategoryBias;
    bool            doGivenBias;
    bool            doGlobalBias;
    Iteration       iteration;
};

ProcLikelihood::Registry registry("ProcLikelihood");

ProcLikelihood::ProcLikelihood(const char *name, const AtomicId *id,
                               MVATrainer *trainer) :
    TrainProcessor(name, id, trainer),
    categoryIdx(-1),
    logOutput(false),
    individual(false),
    neverUndefined(true),
    keepEmpty(false),
    nCategories(1),
    doCategoryBias(false),
    doGivenBias(false),
    doGlobalBias(false),
    iteration(ITER_FILL)
{
}

ProcLikelihood::~ProcLikelihood()
{
}

void ProcLikelihood::configure(DOMElement *elem)
{
    int i = 0;
    bool first = true;
    for(DOMNode *node = elem->getFirstChild();
        node; node = node->getNextSibling()) {
        if (node->getNodeType() != DOMNode::ELEMENT_NODE)
            continue;

        DOMElement *elem = static_cast<DOMElement*>(node);

        XMLSimpleStr nodeName(node->getNodeName());

        if (std::strcmp(nodeName, "general") == 0) {
            if (!first)
                throw cms::Exception("ProcLikelihood")
                    << "Config tag general needs to come "
                       "first." << std::endl;

            if (XMLDocument::hasAttribute(elem, "bias")) {
                double globalBias =
                    XMLDocument::readAttribute<double>(
                                elem, "bias");
                bias.push_back(globalBias);
                doGivenBias = true;
            } else
                doGivenBias = false;

            doCategoryBias = XMLDocument::readAttribute<bool>(
                        elem, "category_bias", false);
            doGlobalBias = XMLDocument::readAttribute<bool>(
                        elem, "global_bias", false);
            logOutput = XMLDocument::readAttribute<bool>(
                        elem, "log", false);
            individual = XMLDocument::readAttribute<bool>(
                        elem, "individual", false);
            neverUndefined = !XMLDocument::readAttribute<bool>(
                        elem, "strict", false);
            keepEmpty = !XMLDocument::readAttribute<bool>(
                        elem, "ignore_empty", true);

            first = false;
            continue;
        }
        first = false;

        if (std::strcmp(nodeName, "bias_table") == 0) {
            if (!bias.empty())
                throw cms::Exception("ProcLikelihood")
                    << "Bias can be only specified once."
                    << std::endl;

            for(DOMNode *subNode = node->getFirstChild();
                subNode; subNode = subNode->getNextSibling()) {
                if (subNode->getNodeType() !=
                    DOMNode::ELEMENT_NODE)
                    continue;

                if (std::strcmp(XMLSimpleStr(
                        subNode->getNodeName()),
                                "bias") != 0)
                    throw cms::Exception("ProcLikelihood")
                        << "Expected bias tag in "
                           "config." << std::endl;

                bias.push_back(
                    XMLDocument::readContent<double>(
                                subNode));
            }

            continue;
        }

        if (std::strcmp(nodeName, "category") != 0) {
            i++;
            continue;
        }

        if (categoryIdx >= 0)
            throw cms::Exception("ProcLikelihood")
                << "More than one category variable given."
                << std::endl;


        unsigned int count = XMLDocument::readAttribute<unsigned int>(
                                elem, "count");

        categoryIdx = i;
        nCategories = count;
    }

    for(DOMNode *node = elem->getFirstChild();
        node; node = node->getNextSibling()) {
        if (node->getNodeType() != DOMNode::ELEMENT_NODE)
            continue;

        XMLSimpleStr nodeName(node->getNodeName());
        if (std::strcmp(nodeName, "general") == 0 ||
            std::strcmp(nodeName, "bias_table") == 0 ||
            std::strcmp(nodeName, "category") == 0)
            continue;

        if (std::strcmp(nodeName, "sigbkg") != 0)
            throw cms::Exception("ProcLikelihood")
                << "Expected sigbkg tag in config section."
                << std::endl;
        elem = static_cast<DOMElement*>(node);

        SigBkg pdf;

        unsigned int size = XMLDocument::readAttribute<unsigned int>(
                            elem, "size", 50);
        pdf.signal.distr.resize(size);
        pdf.background.distr.resize(size);

        pdf.smooth = XMLDocument::readAttribute<unsigned int>(
                            elem, "smooth", 0);

        if (XMLDocument::hasAttribute(elem, "lower") &&
            XMLDocument::hasAttribute(elem, "upper")) {
            pdf.signal.range.min =
                XMLDocument::readAttribute<double>(
                                elem, "lower");
            pdf.signal.range.max =
                XMLDocument::readAttribute<double>(
                                elem, "upper");
            pdf.background.range = pdf.signal.range;
            pdf.iteration = ITER_FILL;
        } else
            pdf.iteration = ITER_EMPTY;

        for(unsigned int i = 0; i < nCategories; i++)
            pdfs.push_back(pdf);
    }

    unsigned int nInputs = getInputs().size();
    if (categoryIdx >= 0)
        nInputs--;

    sigSum.resize(nCategories);
    bkgSum.resize(nCategories);

    if (!doGivenBias && !bias.empty()) {
        doGivenBias = true;
        if (bias.size() != nCategories)
            throw cms::Exception("ProcLikelihood")
                << "Invalid number of category bias entries."
                << std::endl;
    }
    while (doGivenBias && bias.size() < nCategories)
        bias.push_back(bias.front());

    if (pdfs.size() != nInputs * nCategories)
        throw cms::Exception("ProcLikelihood")
            << "Got " << (pdfs.size() / nCategories)
                << " pdf configs for " << nInputs
                << " input variables." << std::endl;
}

Calibration::VarProcessor *ProcLikelihood::getCalibration() const
{
    typedef Calibration::ProcLikelihood Calib;

    Calibration::ProcLikelihood *calib = new Calibration::ProcLikelihood;

    std::vector<unsigned int> pdfMap;
    for(unsigned int i = 0; i < nCategories; i++)
        for(unsigned int j = i; j < pdfs.size(); j += nCategories)
            pdfMap.push_back(j);

    double totalSig = std::accumulate(sigSum.begin(), sigSum.end(), 0.0);
    double totalBkg = std::accumulate(bkgSum.begin(), bkgSum.end(), 0.0);

    for(unsigned int i = 0; i < pdfs.size(); i++) {
        const SigBkg *iter = &pdfs[pdfMap[i]];
        Calibration::ProcLikelihood::SigBkg pdf;

        pdf.signal = Calibration::HistogramF(iter->signal.distr.size(),
                                             iter->signal.range);
        double factor = std::accumulate(iter->signal.distr.begin(),
                                    iter->signal.distr.end(), 0.0);
        if (factor < 1e-20)
            factor = 1.0;
        else
            factor = 1.0 / factor;
        std::vector<double> values(iter->signal.distr.size() + 2);
        std::transform(iter->signal.distr.begin(),
                       iter->signal.distr.end(),
                       values.begin() + 1,
                       [&](auto const& c) {return c * factor;});
        pdf.signal.setValues(values);

        pdf.background =
            Calibration::HistogramF(iter->background.distr.size(),
                                    iter->background.range.min,
                                    iter->background.range.max);
        factor = std::accumulate(iter->background.distr.begin(),
                                 iter->background.distr.end(), 0.0);
        if (factor < 1e-20)
            factor = 1.0;
        else
            factor = 1.0 / factor;
        std::transform(iter->background.distr.begin(),
                       iter->background.distr.end(),
                       values.begin() + 1,
                       [&](auto const& c){return c * factor;});
        pdf.background.setValues(values);

        pdf.useSplines = true;

        calib->pdfs.push_back(pdf);
    }

    calib->categoryIdx = categoryIdx;
    calib->logOutput = logOutput;
    calib->individual = individual;
    calib->neverUndefined = neverUndefined;
    calib->keepEmpty = keepEmpty;

    if (doGlobalBias || doCategoryBias || doGivenBias) {
        for(unsigned int i = 0; i < nCategories; i++) {
            double bias = doGlobalBias
                        ? totalSig / totalBkg
                        : 1.0;
            if (doGivenBias)
                bias *= this->bias[i];
            if (doCategoryBias)
                bias *= (sigSum[i] / totalSig) /
                        (bkgSum[i] / totalBkg);
            calib->bias.push_back(bias);
        }
    }

    return calib;
}

void ProcLikelihood::trainBegin()
{
}

void ProcLikelihood::trainData(const std::vector<double> *values,
                               bool target, double weight)
{
    int category = 0;
    if (categoryIdx >= 0)
        category = (int)values[categoryIdx].front();
    if (category < 0 || category >= (int)nCategories)
        return;

    if (iteration == ITER_FILL) {
        if (target)
            sigSum[category] += weight;
        else
            bkgSum[category] += weight;
    }

    int i = 0;
    for(std::vector<SigBkg>::iterator iter = pdfs.begin() + category;
        iter < pdfs.end(); iter += nCategories, values++) {
        if (i++ == categoryIdx)
            values++;

        switch(iter->iteration) {
            case ITER_EMPTY:
            for(std::vector<double>::const_iterator value =
                            values->begin();
                value != values->end(); value++) {
                iter->signal.range.min =
                    iter->signal.range.max = *value;
                iter->iteration = ITER_RANGE;
                break;
            }
            case ITER_RANGE:
            for(std::vector<double>::const_iterator value =
                            values->begin();
                value != values->end(); value++) {
                iter->signal.range.min =
                    std::min(iter->signal.range.min,
                             *value);
                iter->signal.range.max =
                    std::max(iter->signal.range.max,
                             *value);
            }
            continue;
            case ITER_FILL:
            break;
            default:
            continue;
        }

        PDF &pdf = target ? iter->signal : iter->background;
        unsigned int n = pdf.distr.size() - 1;
        double mult = 1.0 / pdf.range.width();
 
        for(std::vector<double>::const_iterator value =
            values->begin(); value != values->end(); value++) {
            double x = (*value - pdf.range.min) * mult;
            if (x < 0.0)
                x = 0.0;
            else if (x >= 1.0)
                x = 1.0;

            pdf.distr[(unsigned int)(x * n + 0.5)] += weight;
        }
    }
}

void smoothArray(unsigned int n, double *values, unsigned int nTimes)
{
    for(unsigned int iter = 0; iter < nTimes; iter++) {
        double hold = n > 0 ? values[0] : 0.0;
        for(unsigned int i = 0; i < n; i++) {
            double delta = hold * 0.1;
            double rem = 0.0;
            if (i > 0) {
                values[i - 1] += delta;
                rem -= delta;
            }
            if (i < n - 1) {
                hold = values[i + 1];
                values[i + 1] += delta;
                rem -= delta;
            }
            values[i] += rem;
        }
    }
}

void ProcLikelihood::trainEnd()
{
    bool done = true;
    if (iteration == ITER_FILL)
        iteration = ITER_DONE;

    for(std::vector<SigBkg>::iterator iter = pdfs.begin();
        iter != pdfs.end(); iter++) {
        switch(iter->iteration) {
            case ITER_EMPTY:
            case ITER_RANGE:
            iter->background.range = iter->signal.range;
            iter->iteration = ITER_FILL;
            done = false;
            break;
            case ITER_FILL:
            iter->signal.distr.front() *= 2;
            iter->signal.distr.back() *= 2;
            smoothArray(iter->signal.distr.size(),
                        &iter->signal.distr.front(),
                        iter->smooth);

            iter->background.distr.front() *= 2;
            iter->background.distr.back() *= 2;
            smoothArray(iter->background.distr.size(),
                        &iter->background.distr.front(),
                        iter->smooth);

            iter->iteration = ITER_DONE;
            break;
            default:
            /* shut up */;
        }
    }

    if (done)
        trained = true;

    if (done && monitoring) {
        std::vector<SourceVariable*> inputs = getInputs().get();
        if (categoryIdx >= 0)
            inputs.erase(inputs.begin() + categoryIdx);

        for(std::vector<SigBkg>::iterator iter = pdfs.begin();
            iter != pdfs.end(); iter++) {
            unsigned int idx = iter - pdfs.begin();
            unsigned int catIdx = idx % nCategories;
            unsigned int varIdx = idx / nCategories;
            SourceVariable *var = inputs[varIdx];
            std::string name =
                (const char*)var->getSource()->getName()
                + std::string("_")
                + (const char*)var->getName();
            std::string title = name;
            if (categoryIdx >= 0) {
                name += Form("_CAT%d", catIdx);
                title += Form(" (cat. %d)", catIdx);
            }

            unsigned int n = iter->signal.distr.size() - 1;
            double min = iter->signal.range.min -
                         0.5 * iter->signal.range.width() / n;
            double max = iter->signal.range.max +
                         0.5 * iter->signal.range.width() / n;
            TH1F *histo = monitoring->book<TH1F>(name + "_sig",
                (name + "_sig").c_str(),
                (title + " signal").c_str(), n + 1, min, max);
            for(unsigned int i = 0; i < n; i++)
                histo->SetBinContent(
                    i + 1, iter->signal.distr[i]);

            n = iter->background.distr.size() - 1;
            min = iter->background.range.min -
                  0.5 * iter->background.range.width() / n;
            max = iter->background.range.max +
                  0.5 * iter->background.range.width() / n;
            histo = monitoring->book<TH1F>(name + "_bkg",
                (name + "_bkg").c_str(),
                (title + " background").c_str(),
                n + 1, min, max);
            for(unsigned int i = 0; i < n; i++)
                histo->SetBinContent(
                    i + 1, iter->background.distr[i]);
        }
    }
}

void xmlParsePDF(ProcLikelihood::PDF &pdf, DOMElement *elem)
{
    if (!elem ||
        std::strcmp(XMLSimpleStr(elem->getNodeName()), "pdf") != 0)
        throw cms::Exception("ProcLikelihood")
            << "Expected pdf tag in sigbkg train data."
            << std::endl;

    pdf.range.min = XMLDocument::readAttribute<double>(elem, "lower");
    pdf.range.max = XMLDocument::readAttribute<double>(elem, "upper");

    pdf.distr.clear();
    for(DOMNode *node = elem->getFirstChild();
        node; node = node->getNextSibling()) {
        if (node->getNodeType() != DOMNode::ELEMENT_NODE)
            continue;

        if (std::strcmp(XMLSimpleStr(node->getNodeName()),
                                     "value") != 0)
            throw cms::Exception("ProcLikelihood")
                << "Expected value tag in train file."
                << std::endl;

        pdf.distr.push_back(XMLDocument::readContent<double>(node));
    }
}

namespace {
    struct Id {
        AtomicId    source;
        AtomicId    name;
        unsigned int    category;

        inline Id(AtomicId source, AtomicId name,
                  unsigned int category) :
            source(source), name(name), category(category) {}

        inline bool operator == (const Id &other) const
        {
            return source == other.source &&
                   name == other.name &&
                   category == other.category;
        }

        inline bool operator < (const Id &other) const
        {
            if (source < other.source)
                return true;
            if (!(source == other.source))
                return false;
            if (name < other.name)
                return true;
            if (!(name == other.name))
                return false;
            return category < other.category;
        }
    };
}

bool ProcLikelihood::load()
{
    std::string filename = trainer->trainFileName(this, "xml");
    if (!exists(filename))
        return false;

    XMLDocument xml(filename);
    DOMElement *elem = xml.getRootNode();
    if (std::strcmp(XMLSimpleStr(elem->getNodeName()),
                                 "ProcLikelihood") != 0)
        throw cms::Exception("ProcLikelihood")
            << "XML training data file has bad root node."
            << std::endl;

    unsigned int version = XMLDocument::readAttribute<unsigned int>(
                            elem, "version", 1);

    if (version < 1 || version > 2)
        throw cms::Exception("ProcLikelihood")
            << "Unsupported version " << version
            << "in train file." << std::endl;

    DOMNode *node;
    for(node = elem->getFirstChild();
        node; node = node->getNextSibling()) {
        if (node->getNodeType() != DOMNode::ELEMENT_NODE)
            continue;

        if (std::strcmp(XMLSimpleStr(node->getNodeName()),
                        "categories") != 0)
            throw cms::Exception("ProcLikelihood")
                << "Expected categories tag in train file."
                << std::endl;

        unsigned int i = 0;
        for(DOMNode *subNode = node->getFirstChild();
            subNode; subNode = subNode->getNextSibling()) {
            if (subNode->getNodeType() != DOMNode::ELEMENT_NODE)
                continue;

            if (i >= nCategories)
                throw cms::Exception("ProcLikelihood")
                    << "Too many categories in train "
                           "file." << std::endl;

            if (std::strcmp(XMLSimpleStr(subNode->getNodeName()),
                            "category") != 0)
                throw cms::Exception("ProcLikelihood")
                    << "Expected category tag in train "
                           "file." << std::endl;

            elem = static_cast<DOMElement*>(subNode);

            sigSum[i] = XMLDocument::readAttribute<double>(
                            elem, "signal");
            bkgSum[i] = XMLDocument::readAttribute<double>(
                            elem, "background");
            i++;
        }
        if (i < nCategories)
            throw cms::Exception("ProcLikelihood")
                << "Too few categories in train file."
                << std::endl;

        break;
    }

    std::map<Id, SigBkg*> pdfMap;

    for(std::vector<SigBkg>::iterator iter = pdfs.begin();
        iter != pdfs.end(); ++iter) {
        SigBkg *ptr = &*iter;
        unsigned int idx = iter - pdfs.begin();
        unsigned int catIdx = idx % nCategories;
        unsigned int varIdx = idx / nCategories;
        if (categoryIdx >= 0 && (int)varIdx >= categoryIdx)
            varIdx++;
        const SourceVariable *var = getInputs().get()[varIdx];
        Id id(var->getSource()->getName(), var->getName(), catIdx);

        pdfMap[id] = ptr;
    }

    std::vector<SigBkg>::iterator cur = pdfs.begin();

    for(node = node->getNextSibling();
        node; node = node->getNextSibling()) {
        if (node->getNodeType() != DOMNode::ELEMENT_NODE)
            continue;

        if (std::strcmp(XMLSimpleStr(node->getNodeName()),
                        "sigbkg") != 0)
            throw cms::Exception("ProcLikelihood")
                << "Expected sigbkg tag in train file."
                << std::endl;
        elem = static_cast<DOMElement*>(node);

        SigBkg *pdf = nullptr;
        switch(version) {
            case 1:
            if (cur == pdfs.end())
                throw cms::Exception("ProcLikelihood")
                    << "Superfluous SigBkg in train data."
                    << std::endl;
            pdf = &*cur++;
            break;
            case 2: {
            Id id(XMLDocument::readAttribute<std::string>(
                            elem, "source"),
                  XMLDocument::readAttribute<std::string>(
                            elem, "name"),
                  XMLDocument::readAttribute<unsigned int>(
                                                        elem, "category", 0));
            std::map<Id, SigBkg*>::const_iterator pos =
                            pdfMap.find(id);
            if (pos == pdfMap.end())
                continue;
            else
                pdf = pos->second;
            }   break;
        }

        for(node = elem->getFirstChild();
            node && node->getNodeType() != DOMNode::ELEMENT_NODE;
            node = node->getNextSibling());
        DOMElement *elemSig =
                node ? static_cast<DOMElement*>(node) : nullptr;

        for(node = node->getNextSibling();
            node && node->getNodeType() != DOMNode::ELEMENT_NODE;
            node = node->getNextSibling());
        while(node && node->getNodeType() != DOMNode::ELEMENT_NODE)
            node = node->getNextSibling();
        DOMElement *elemBkg =
                node ? static_cast<DOMElement*>(node) : nullptr;

        for(node = node->getNextSibling();
            node && node->getNodeType() != DOMNode::ELEMENT_NODE;
            node = node->getNextSibling());
        if (node)
            throw cms::Exception("ProcLikelihood")
                << "Superfluous tags in sigbkg train data."
                << std::endl;

        xmlParsePDF(pdf->signal, elemSig);
        xmlParsePDF(pdf->background, elemBkg);

        pdf->iteration = ITER_DONE;

        node = elem;
    }

    if (version == 1 && cur != pdfs.end())
        throw cms::Exception("ProcLikelihood")
            << "Missing SigBkg in train data." << std::endl;

    iteration = ITER_DONE;
    trained = true;
    for(std::vector<SigBkg>::const_iterator iter = pdfs.begin();
        iter != pdfs.end(); ++iter) {
        if (iter->iteration != ITER_DONE) {
            trained = false;
            break;
        }
    }

    return true;
}

DOMElement *xmlStorePDF(DOMDocument *doc,
                               const ProcLikelihood::PDF &pdf)
{
    DOMElement *elem = doc->createElement(XMLUniStr("pdf"));

    XMLDocument::writeAttribute(elem, "lower", pdf.range.min);
    XMLDocument::writeAttribute(elem, "upper", pdf.range.max);

    for(std::vector<double>::const_iterator iter =
        pdf.distr.begin(); iter != pdf.distr.end(); iter++) {
        DOMElement *value = doc->createElement(XMLUniStr("value"));
        elem->appendChild(value);   

        XMLDocument::writeContent<double>(value, doc, *iter);
    }

    return elem;
}

void ProcLikelihood::save()
{
    XMLDocument xml(trainer->trainFileName(this, "xml"), true);
    DOMDocument *doc = xml.createDocument("ProcLikelihood");
    XMLDocument::writeAttribute(doc->getDocumentElement(), "version", 2);

    DOMElement *elem = doc->createElement(XMLUniStr("categories"));
    xml.getRootNode()->appendChild(elem);
    for(unsigned int i = 0; i < nCategories; i++) {
        DOMElement *category = doc->createElement(XMLUniStr("category"));
        elem->appendChild(category);
        XMLDocument::writeAttribute(category, "signal", sigSum[i]);
        XMLDocument::writeAttribute(category, "background", bkgSum[i]);
    }

    for(std::vector<SigBkg>::const_iterator iter = pdfs.begin();
        iter != pdfs.end(); iter++) {
        elem = doc->createElement(XMLUniStr("sigbkg"));
        xml.getRootNode()->appendChild(elem);

        unsigned int idx = iter - pdfs.begin();
        unsigned int catIdx = idx % nCategories;
        unsigned int varIdx = idx / nCategories;
        if (categoryIdx >= 0 && (int)varIdx >= categoryIdx)
            varIdx++;
        const SourceVariable *var = getInputs().get()[varIdx];
        XMLDocument::writeAttribute(elem, "source",
                (const char*)var->getSource()->getName());
        XMLDocument::writeAttribute(elem, "name",
                (const char*)var->getName());
        if (categoryIdx >= 0)
            XMLDocument::writeAttribute(elem, "category", catIdx);

        elem->appendChild(xmlStorePDF(doc, iter->signal));
        elem->appendChild(xmlStorePDF(doc, iter->background));
    }
}

} // anonymous namespace