VarProcessor.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:     MVAComputer
// Class  :     VarProcessor
// 

// Implementation:
//     Base class for variable processors. Basically only passes calls
//     through to virtual methods in the actual implementation daughter class.
//
// Author:      Christophe Saout
// Created:     Sat Apr 24 15:18 CEST 2007
// $Id: VarProcessor.cc,v 1.7 2010/02/19 10:28:37 saout Exp $
//

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

#include "FWCore/PluginManager/interface/PluginManager.h"
#include "FWCore/PluginManager/interface/PluginFactory.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "PhysicsTools/MVAComputer/interface/VarProcessor.h"
#include "PhysicsTools/MVAComputer/interface/Calibration.h"
#include "PhysicsTools/MVAComputer/interface/BitSet.h"

// #define DEBUG_DERIV

#ifdef DEBUG_DERIV
#   include <Reflex/Tools.h>
#endif

EDM_REGISTER_PLUGINFACTORY(PhysicsTools::VarProcessor::PluginFactory,
                           "PhysicsToolsMVAComputer");

namespace PhysicsTools {

VarProcessor::VarProcessor(const char *name,
                           const Calibration::VarProcessor *calib,
                           const MVAComputer *computer) :
    computer(computer),
    inputVars(Calibration::convert(calib->inputVars)),
    nInputVars(inputVars.bits())
{
}

VarProcessor::~VarProcessor()
{
    inputVars = BitSet(0);
    nInputVars = 0;
}

void VarProcessor::configure(ConfigCtx &config)
{
    ConfigCtx::size_type pos = config.size();
    if (pos != inputVars.size())
        return;

    ConfIterator iter(inputVars.iter(), config);
    configure(iter, nInputVars);

    VarProcessor *loop = config.loop ? config.loop : this;
    ConfigCtx::Context *ctx =
        loop->configureLoop(config.ctx, config.begin(),
                            config.begin() + pos, config.end());

    if (ctx != config.ctx) {
        delete config.ctx;
        config.ctx = ctx;
    }

    if (config.loop && !ctx)
        config.loop = 0;
    else if (!config.loop && ctx)
        config.loop = this;
}

VarProcessor::ConfigCtx::ConfigCtx(std::vector<Variable::Flags> flags) :
    loop(0), ctx(0)
{
    for(std::vector<Variable::Flags>::const_iterator iter = flags.begin();
        iter != flags.end(); ++iter)
        configs.push_back(Config(*iter, 1));
}

VarProcessor::ConfigCtx::Context *
VarProcessor::configureLoop(ConfigCtx::Context *ctx, ConfigCtx::iterator begin,
                            ConfigCtx::iterator cur, ConfigCtx::iterator end)
{
    return 0;
}

template<>
VarProcessor *ProcessRegistry<VarProcessor, Calibration::VarProcessor,
                              const MVAComputer>::Factory::create(
            const char *name, const Calibration::VarProcessor *calib,
        const MVAComputer *parent)
{
    VarProcessor *result = ProcessRegistry::create(name, calib, parent);
    if (!result) {
        // try to load the shared library and retry
        try {
            delete VarProcessor::PluginFactory::get()->create(
                    std::string("VarProcessor/") + name);
            result = ProcessRegistry::create(name, calib, parent);
        } catch(const cms::Exception &e) {
            // caller will have to deal with the null pointer
            // in principle this will just give a slightly more
            // descriptive error message (and will rethrow anyhow)

                  edm::LogError("CannotBuildMVAProc")
                    << "Caught exception when building processor: "
                    << name << " message: " << std::endl
                    << e.what() << std::endl;
                  throw e;
        }
    }
    return result;
}

void VarProcessor::deriv(double *input, int *conf, double *output,
                         int *outConf, int *loop, unsigned int offset,
                         unsigned int in, unsigned int out_,
                         std::vector<double> &deriv) const
{
    ValueIterator iter(inputVars.iter(), input, conf,
                       output, outConf, loop, offset);

    eval(iter, nInputVars);

    std::vector<double> matrix = this->deriv(iter, nInputVars);

    unsigned int size = 0;
    while(iter)
        size += (iter++).size();
    unsigned int out = matrix.empty() ? 0 : (matrix.size() / size);

    if (matrix.size() != out * size ||
        (out > 1 && (int)out != outConf[out_] - outConf[0]))
        throw cms::Exception("VarProcessor")
            << "Derivative matrix implausible size in "
            << typeid(*this).name() << "."
            << std::endl;

#ifdef DEBUG_DERIV
    if (!matrix.empty()) {
        std::cout << "---------------- "
                  << ROOT::Reflex::Tools::Demangle(typeid(*this))
                  << std::endl;
        for(unsigned int i = 0; i < out; i++) {
            for(unsigned int j = 0; j < size; j++)
                std::cout << matrix.at(i*size+j) << "\t";
            std::cout << std::endl;
        }
        std::cout << "----------------" << std::endl;
    }

    std::cout << "======= in = " << in << ", size = " << size
              << ", out = " << out << ", matrix = " << matrix.size()
              << std::endl;
#endif

    unsigned int sz = (outConf[out_] - in) * in;
    unsigned int oldSz = deriv.size();
    if (oldSz < sz)
        deriv.resize(sz);

    double *begin = &deriv.front() + (outConf[0] - in + offset) * in;
    double *end = begin + out * in;
    if (begin < &deriv.front() + oldSz)
        std::fill(begin, end, 0.0);

    if (matrix.empty())
        return;

    double *m0 = &matrix.front();
    BitSet::Iterator cur = inputVars.iter();
    for(unsigned int i = 0; i < nInputVars; i++, ++cur) {
#ifdef DEBUG_DERIV
        std::cout << " inputvar " << i << std::endl;
#endif
        int *curConf = conf + cur();
        unsigned int pos = *curConf;
#ifdef DEBUG_DERIV
        std::cout << " -> cur = " << cur() << ", pos = "
                  << pos << std::endl;
#endif
        if (loop && curConf >= loop) {
            pos += offset;
            loop = 0;
        }

        unsigned int n = loop ? (curConf[1] - curConf[0]) : 1;
        for(unsigned int j = 0; j < n; m0++, j++, pos++) {
#ifdef DEBUG_DERIV
            std::cout << "  multip " << j << std::endl;
#endif
            double *p = begin;
            if (pos >= in) {
#ifdef DEBUG_DERIV
                std::cout << "   deriv " << (pos - in)
                          << std::endl;
#endif
                const double *q = &deriv.front() +
                                  (pos - in) * in;
                for(const double *m = m0; p < end;
                    m += size, q -= in)
                    for(unsigned int k = 0; k < in; k++)
                        *p++ += *m * *q++;
            } else {
#ifdef DEBUG_DERIV
                std::cout << "   in " << pos << std::endl;
#endif
                for(const double *m = m0; p < end;
                    m += size, p += in)
                    p[pos] += *m;
            }
        }
    }

#ifdef DEBUG_DERIV
    std::cout << "================" << std::endl;
    for(const double *p = &deriv.front();
        p != &deriv.front() + deriv.size();) {
        for(unsigned int j = 0; j < in; j++)
            std::cout << *p++ << "\t";
        std::cout << std::endl;
    }
    std::cout << "================" << std::endl;
#endif
}

} // namespace PhysicsTools