MVAComputer/src/VarProcessor.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:     MVAComputer
0004 // Class  :     VarProcessor
0005 //
0006
0007 // Implementation:
0008 //     Base class for variable processors. Basically only passes calls
0009 //     through to virtual methods in the actual implementation daughter class.
0010 //
0011 // Author:      Christophe Saout
0012 // Created:     Sat Apr 24 15:18 CEST 2007
0013 //
0014
0015 #include <cassert>
0016
0017 #include "FWCore/Utilities/interface/Exception.h"
0018
0019 #include "FWCore/PluginManager/interface/PluginManager.h"
0020 #include "FWCore/PluginManager/interface/PluginFactory.h"
0021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0022
0023 #include "PhysicsTools/MVAComputer/interface/VarProcessor.h"
0024 #include "PhysicsTools/MVAComputer/interface/Calibration.h"
0025 #include "PhysicsTools/MVAComputer/interface/BitSet.h"
0026 #include "PhysicsTools/MVAComputer/interface/ProcessRegistry.icc"
0027
0028 // #define DEBUG_DERIV
0029
0030 #ifdef DEBUG_DERIV
0031 #include "FWCore/Utilities/interface/TypeDemangler.h"
0032 #endif
0033
0034 typedef edmplugin::PluginFactory<PhysicsTools::VarProcessor::PluginFunctionPrototype> VPPluginFactory;
0035 EDM_REGISTER_PLUGINFACTORY(VPPluginFactory, "PhysicsToolsMVAComputer");
0036
0037 namespace PhysicsTools {
0038
0039   VarProcessor::VarProcessor(const char *name, const Calibration::VarProcessor *calib, const MVAComputer *computer)
0040       : computer(computer), inputVars(Calibration::convert(calib->inputVars)), nInputVars(inputVars.bits()) {}
0041
0042   VarProcessor::~VarProcessor() {
0043     inputVars = BitSet(0);
0044     nInputVars = 0;
0045   }
0046
0047   void VarProcessor::configure(ConfigCtx &config) {
0048     ConfigCtx::size_type pos = config.size();
0049     if (pos != inputVars.size())
0050       return;
0051
0052     ConfIterator iter(inputVars.iter(), config);
0053     configure(iter, nInputVars);
0054
0055     VarProcessor *loop = config.loop ? config.loop : this;
0056     ConfigCtx::Context *ctx = loop->configureLoop(config.ctx, config.begin(), config.begin() + pos, config.end());
0057
0058     if (ctx != config.ctx) {
0059       delete config.ctx;
0060       config.ctx = ctx;
0061     }
0062
0063     if (config.loop && !ctx)
0064       config.loop = nullptr;
0065     else if (!config.loop && ctx)
0066       config.loop = this;
0067   }
0068
0069   VarProcessor::ConfigCtx::ConfigCtx(const std::vector<Variable::Flags> &flags) : loop(nullptr), ctx(nullptr) {
0070     for (std::vector<Variable::Flags>::const_iterator iter = flags.begin(); iter != flags.end(); ++iter)
0071       configs.push_back(Config(*iter, 1));
0072   }
0073
0074   VarProcessor::ConfigCtx::Context *VarProcessor::configureLoop(ConfigCtx::Context *ctx,
0075                                                                 ConfigCtx::iterator begin,
0076                                                                 ConfigCtx::iterator cur,
0077                                                                 ConfigCtx::iterator end) {
0078     return nullptr;
0079   }
0080
0081   template <>
0082   VarProcessor *ProcessRegistry<VarProcessor, Calibration::VarProcessor, const MVAComputer>::Factory::create(
0083       const char *name, const Calibration::VarProcessor *calib, const MVAComputer *parent) {
0084     VarProcessor *result = ProcessRegistry::create(name, calib, parent);
0085     if (!result) {
0086       // try to load the shared library and retry
0087       try {
0088         std::unique_ptr<PhysicsTools::VarProcessor::Dummy> tmp{
0089             VPPluginFactory::get()->create(std::string("VarProcessor/") + name)};
0090         result = ProcessRegistry::create(name, calib, parent);
0091       } catch (const cms::Exception &e) {
0092         // caller will have to deal with the null pointer
0093         // in principle this will just give a slightly more
0094         // descriptive error message (and will rethrow anyhow)
0095
0096         edm::LogError("CannotBuildMVAProc")
0097             << "Caught exception when building processor: " << name << " message: " << std::endl
0098             << e.what() << std::endl;
0099         throw e;
0100       }
0101     }
0102     return result;
0103   }
0104
0105   void VarProcessor::deriv(double *input,
0106                            int *conf,
0107                            double *output,
0108                            int *outConf,
0109                            int *loop,
0110                            LoopCtx &ctx,
0111                            unsigned int offset,
0112                            unsigned int in,
0113                            unsigned int out_,
0114                            std::vector<double> &deriv) const {
0115     ValueIterator iter(inputVars.iter(), input, conf, output, outConf, loop, ctx, offset);
0116
0117     eval(iter, nInputVars);
0118
0119     std::vector<double> matrix = this->deriv(iter, nInputVars);
0120
0121     unsigned int size = 0;
0122     while (iter)
0123       size += (iter++).size();
0124     bool empty = matrix.empty();
0125     assert(size != 0 || empty);
0126     unsigned int out = empty ? 0 : (matrix.size() / size);
0127
0128     if (matrix.size() != out * size || (out > 1 && (int)out != outConf[out_] - outConf[0]))
0129       throw cms::Exception("VarProcessor")
0130           << "Derivative matrix implausible size in " << typeid(*this).name() << "." << std::endl;
0131
0132 #ifdef DEBUG_DERIV
0133     if (!matrix.empty()) {
0134       std::string demangledName;
0135       edm::typeDemangle(typeid(*this).name(), demangledName);
0136       std::cout << demangledName << std::endl;
0137       for (unsigned int i = 0; i < out; i++) {
0138         for (unsigned int j = 0; j < size; j++)
0139           std::cout << matrix.at(i * size + j) << "\t";
0140         std::cout << std::endl;
0141       }
0142       std::cout << "----------------" << std::endl;
0143     }
0144
0145     std::cout << "======= in = " << in << ", size = " << size << ", out = " << out << ", matrix = " << matrix.size()
0146               << std::endl;
0147 #endif
0148
0149     unsigned int sz = (outConf[out_] - in) * in;
0150     unsigned int oldSz = deriv.size();
0151     if (oldSz < sz)
0152       deriv.resize(sz);
0153
0154     double *begin = &deriv.front() + (outConf[0] - in + offset) * in;
0155     double *end = begin + out * in;
0156     if (begin < &deriv.front() + oldSz)
0157       std::fill(begin, end, 0.0);
0158
0159     if (matrix.empty())
0160       return;
0161
0162     double *m0 = &matrix.front();
0163     BitSet::Iterator cur = inputVars.iter();
0164     for (unsigned int i = 0; i < nInputVars; i++, ++cur) {
0165 #ifdef DEBUG_DERIV
0166       std::cout << " inputvar " << i << std::endl;
0167 #endif
0168       int *curConf = conf + cur();
0169       unsigned int pos = *curConf;
0170 #ifdef DEBUG_DERIV
0171       std::cout << " -> cur = " << cur() << ", pos = " << pos << std::endl;
0172 #endif
0173       if (loop && curConf >= loop) {
0174         pos += offset;
0175         loop = nullptr;
0176       }
0177
0178       unsigned int n = loop ? (curConf[1] - curConf[0]) : 1;
0179       for (unsigned int j = 0; j < n; m0++, j++, pos++) {
0180 #ifdef DEBUG_DERIV
0181         std::cout << "  multip " << j << std::endl;
0182 #endif
0183         double *p = begin;
0184         if (pos >= in) {
0185 #ifdef DEBUG_DERIV
0186           std::cout << "   deriv " << (pos - in) << std::endl;
0187 #endif
0188           const double *q = &deriv.front() + (pos - in) * in;
0189           for (const double *m = m0; p < end; m += size, q -= in)
0190             for (unsigned int k = 0; k < in; k++)
0191               *p++ += *m * *q++;
0192         } else {
0193 #ifdef DEBUG_DERIV
0194           std::cout << "   in " << pos << std::endl;
0195 #endif
0196           for (const double *m = m0; p < end; m += size, p += in)
0197             p[pos] += *m;
0198         }
0199       }
0200     }
0201
0202 #ifdef DEBUG_DERIV
0203     std::cout << "================" << std::endl;
0204     for (const double *p = &deriv.front(); p != &deriv.front() + deriv.size();) {
0205       for (unsigned int j = 0; j < in; j++)
0206         std::cout << *p++ << "\t";
0207       std::cout << std::endl;
0208     }
0209     std::cout << "================" << std::endl;
0210 #endif
0211   }
0212
0213 }  // namespace PhysicsTools
0214
0215 // Force instantiation of its templated static methods.
0216 template void PhysicsTools::ProcessRegistry<PhysicsTools::VarProcessor,
0217                                             PhysicsTools::Calibration::VarProcessor,
0218                                             PhysicsTools::MVAComputer const>::unregisterProcess(char const *);
0219 template void PhysicsTools::ProcessRegistry<PhysicsTools::VarProcessor,
0220                                             PhysicsTools::Calibration::VarProcessor,
0221                                             PhysicsTools::MVAComputer const>::
0222     registerProcess(char const *,
0223                     PhysicsTools::ProcessRegistry<PhysicsTools::VarProcessor,
0224                                                   PhysicsTools::Calibration::VarProcessor,
0225                                                   PhysicsTools::MVAComputer const> const *);