Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-11-15 04:15:14

0001 #include <iostream>
0002 #include <sstream>
0003 #include <fstream>
0004 #include <algorithm>
0005 #include <iterator>
0006 #include <cstddef>
0007 #include <cstring>
0008 #include <vector>
0009 #include <set>
0010 
0011 // ROOT version magic to support TMVA interface changes in newer ROOT
0012 #include <RVersion.h>
0013 
0014 #include <TBufferFile.h>
0015 #include <TClass.h>
0016 
0017 #include "FWCore/Utilities/interface/Exception.h"
0018 
0019 #include "PhysicsTools/MVAComputer/interface/zstream.h"
0020 
0021 #include "PhysicsTools/MVAComputer/interface/MVAComputer.h"
0022 #include "PhysicsTools/MVAComputer/interface/VarProcessor.h"
0023 #include "PhysicsTools/MVAComputer/interface/Calibration.h"
0024 #include "PhysicsTools/MVAComputer/interface/Variable.h"
0025 #include "PhysicsTools/MVAComputer/interface/AtomicId.h"
0026 
0027 // #define DEBUG_EVAL
0028 
0029 #ifdef DEBUG_EVAL
0030 #include "FWCore/Utilities/interface/TypeDemangler.h"
0031 #endif
0032 
0033 #define STANDALONE_HEADER "MVAComputer calibration\n"
0034 
0035 namespace PhysicsTools {
0036 
0037   MVAComputer::MVAComputer(const Calibration::MVAComputer *calib) : nVars(0), output(0) { setup(calib); }
0038 
0039   MVAComputer::MVAComputer(Calibration::MVAComputer *calib, bool owned) : nVars(0), output(0) {
0040     if (owned)
0041       this->owned.reset(calib);
0042     setup(calib);
0043   }
0044 
0045   MVAComputer::MVAComputer(const char *filename) : nVars(0), output(0), owned(readCalibration(filename)) {
0046     setup(owned.get());
0047   }
0048 
0049   MVAComputer::MVAComputer(std::istream &is) : nVars(0), output(0), owned(readCalibration(is)) { setup(owned.get()); }
0050 
0051   void MVAComputer::setup(const Calibration::MVAComputer *calib) {
0052     nVars = calib->inputSet.size();
0053     output = calib->output;
0054 
0055     std::vector<Variable::Flags> flags(nVars, Variable::FLAG_ALL);
0056     const TrainMVAComputerCalibration *trainCalib = dynamic_cast<const TrainMVAComputerCalibration *>(calib);
0057     if (trainCalib)
0058       trainCalib->initFlags(flags);
0059 
0060     VarProcessor::ConfigCtx config(flags);
0061     std::vector<const Calibration::VarProcessor *> processors = calib->getProcessors();
0062 
0063     for (const auto *calibProc : processors) {
0064       std::string name = calibProc->getInstanceName();
0065       VarProcessor *processor = VarProcessor::create(name.c_str(), calibProc, this);
0066       if (!processor)
0067         throw cms::Exception("UnknownProcessor") << name << " could not be instantiated." << std::endl;
0068 
0069       VarProcessor::ConfigCtx::iterator::difference_type pos = config.end() - config.begin();
0070       processor->configure(config);
0071       unsigned int nOutput = (config.end() - config.begin()) - pos;
0072       if (!nOutput)
0073         throw cms::Exception("InvalidProcessor") << name
0074                                                  << " rejected input variable "
0075                                                     "configuration"
0076                                                  << std::endl;
0077 
0078       varProcessors.push_back(Processor(processor, nOutput));
0079     }
0080 
0081     for (VarProcessor::ConfigCtx::iterator iter = config.begin() + nVars; iter != config.end(); iter++) {
0082       VarProcessor::Config *origin = &config[iter->origin];
0083       if (iter->origin >= nVars)
0084         iter->origin = origin->origin;
0085 
0086       if (iter->mask & Variable::FLAG_MULTIPLE) {
0087         iter->mask = (Variable::Flags)(iter->mask & origin->mask);
0088         config[iter->origin].origin++;
0089       }
0090     }
0091 
0092     nVars = config.size();
0093 
0094     if (output >= nVars)
0095       // FIXME || config[output].mask != Variable::FLAG_NONE)
0096       throw cms::Exception("InvalidOutput") << "Output variable at index " << output << " invalid." << std::endl;
0097 
0098     std::set<InputVar> variables;
0099     unsigned int i = 0;
0100     for (std::vector<Calibration::Variable>::const_iterator iter = calib->inputSet.begin();
0101          iter != calib->inputSet.end();
0102          ++iter, i++) {
0103       InputVar var;
0104       var.var = Variable(iter->name, config[i].mask);
0105       var.index = i;
0106       var.multiplicity = 0;
0107       variables.insert(var);
0108     }
0109 
0110     inputVariables.resize(i);
0111 
0112     std::copy(variables.begin(), variables.end(), inputVariables.begin());
0113 
0114     for (unsigned int j = 0; j < i; j++)
0115       inputVariables[j].multiplicity = config[j].origin;
0116   }
0117 
0118   MVAComputer::~MVAComputer() {}
0119 
0120   int MVAComputer::getVariableId(AtomicId name) const {
0121     std::vector<InputVar>::const_iterator pos = std::lower_bound(inputVariables.begin(), inputVariables.end(), name);
0122 
0123     if (pos == inputVariables.end() || pos->var.getName() != name)
0124       return -1;
0125     else
0126       return pos->index;
0127   }
0128 
0129   template <class T>
0130   void MVAComputer::evalInternal(T &ctx) const {
0131     double *output = ctx.values() + ctx.n();
0132     int *outConf = ctx.conf() + inputVariables.size();
0133 
0134 #ifdef DEBUG_EVAL
0135     std::cout << "Input" << std::endl;
0136     double *v = ctx.values();
0137     for (int *o = ctx.conf(); o < outConf; o++) {
0138       std::cout << "\tVar " << (o - ctx.conf()) << std::endl;
0139       for (int i = o[0]; i < o[1]; i++)
0140         std::cout << "\t\t" << *v++ << std::endl;
0141     }
0142 #endif
0143     std::vector<Processor>::const_iterator iter = varProcessors.begin();
0144     while (iter != varProcessors.end()) {
0145       std::vector<Processor>::const_iterator loop = iter;
0146       int *loopOutConf = outConf;
0147       int *loopStart = nullptr;
0148       double *loopOutput = output;
0149       VarProcessor::LoopCtx loopCtx;
0150 
0151       VarProcessor::LoopStatus status = VarProcessor::kNext;
0152       unsigned int offset = 0;
0153       while (status != VarProcessor::kStop) {
0154         std::vector<Processor>::const_iterator next = iter + 1;
0155         unsigned int nextOutput = (next != varProcessors.end()) ? next->nOutput : 0;
0156 
0157 #ifdef DEBUG_EVAL
0158         std::string demangledName;
0159         edm::typeDemangle(typeid(*iter->processor).name(), demangledName);
0160         std::cout << demangledName << std::endl;
0161 #endif
0162         if (status != VarProcessor::kSkip)
0163           ctx.eval(
0164               &*iter->processor, outConf, output, loopStart ? loopStart : loopOutConf, loopCtx, offset, iter->nOutput);
0165 
0166 #ifdef DEBUG_EVAL
0167         for (unsigned int i = 0; i < iter->nOutput; i++, outConf++) {
0168           std::cout << "\tVar " << (outConf - ctx.conf()) << std::endl;
0169           for (int j = outConf[0]; j < outConf[1]; j++)
0170             std::cout << "\t\t" << *output++ << std::endl;
0171         }
0172 #else
0173         int orig = *outConf;
0174         outConf += iter->nOutput;
0175         output += *outConf - orig;
0176 #endif
0177 
0178         status = loop->processor->loop(output, outConf, nextOutput, loopCtx, offset);
0179 
0180         if (status == VarProcessor::kReset) {
0181           outConf = loopOutConf;
0182           output = loopOutput;
0183           loopStart = nullptr;
0184           offset = 0;
0185           iter = loop;
0186         } else {
0187           if (loop == iter)
0188             loopStart = outConf;
0189           iter = next;
0190         }
0191       }
0192     }
0193   }
0194 
0195   Calibration::MVAComputer *MVAComputer::readCalibration(const char *filename) {
0196     std::ifstream file(filename);
0197     return readCalibration(file);
0198   }
0199 
0200   Calibration::MVAComputer *MVAComputer::readCalibration(std::istream &is) {
0201     if (!is.good())
0202       throw cms::Exception("InvalidFileState") << "Stream passed to MVAComputer::readCalibration "
0203                                                   "has an invalid state."
0204                                                << std::endl;
0205 
0206     char header[sizeof STANDALONE_HEADER - 1] = {
0207         0,
0208     };
0209     if (is.readsome(header, sizeof header) != sizeof header ||
0210         std::memcmp(header, STANDALONE_HEADER, sizeof header) != 0)
0211       throw cms::Exception("InvalidFileFormat") << "Stream passed to MVAComputer::readCalibration "
0212                                                    "is not a valid calibration file."
0213                                                 << std::endl;
0214 
0215     TClass *rootClass = TClass::GetClass("PhysicsTools::Calibration::MVAComputer");
0216     if (!rootClass)
0217       throw cms::Exception("DictionaryMissing") << "CondFormats dictionary for "
0218                                                    "PhysicsTools::Calibration::MVAComputer missing"
0219                                                 << std::endl;
0220 
0221     ext::izstream izs(&is);
0222     std::ostringstream ss;
0223     ss << izs.rdbuf();
0224     std::string buf = ss.str();
0225 
0226     TBufferFile buffer(TBuffer::kRead, buf.size(), const_cast<void *>(static_cast<const void *>(buf.c_str())), kFALSE);
0227     buffer.InitMap();
0228 
0229     std::unique_ptr<Calibration::MVAComputer> calib(new Calibration::MVAComputer());
0230     buffer.StreamObject(static_cast<void *>(calib.get()), rootClass);
0231 
0232     return calib.release();
0233   }
0234 
0235   void MVAComputer::writeCalibration(const char *filename, const Calibration::MVAComputer *calib) {
0236     std::ofstream file(filename);
0237     writeCalibration(file, calib);
0238   }
0239 
0240   void MVAComputer::writeCalibration(std::ostream &os, const Calibration::MVAComputer *calib) {
0241     if (!os.good())
0242       throw cms::Exception("InvalidFileState") << "Stream passed to MVAComputer::writeCalibration "
0243                                                   "has an invalid state."
0244                                                << std::endl;
0245 
0246     os << STANDALONE_HEADER;
0247 
0248     TClass *rootClass = TClass::GetClass("PhysicsTools::Calibration::MVAComputer");
0249     if (!rootClass)
0250       throw cms::Exception("DictionaryMissing") << "CondFormats dictionary for "
0251                                                    "PhysicsTools::Calibration::MVAComputer missing"
0252                                                 << std::endl;
0253 
0254     TBufferFile buffer(TBuffer::kWrite);
0255     buffer.StreamObject(const_cast<void *>(static_cast<const void *>(calib)), rootClass);
0256 
0257     ext::ozstream ozs(&os);
0258     ozs.write(buffer.Buffer(), buffer.Length());
0259     ozs.flush();
0260   }
0261 
0262   // instantiate use cases fo MVAComputer::evalInternal
0263 
0264   void MVAComputer::DerivContext::eval(const VarProcessor *proc,
0265                                        int *outConf,
0266                                        double *output,
0267                                        int *loop,
0268                                        VarProcessor::LoopCtx &ctx,
0269                                        unsigned int offset,
0270                                        unsigned int out) const {
0271     proc->deriv(values(), conf(), output, outConf, loop, ctx, offset, n(), out, deriv_);
0272   }
0273 
0274   double MVAComputer::DerivContext::output(unsigned int output, std::vector<double> &derivs) const {
0275     derivs.clear();
0276     derivs.reserve(n_);
0277 
0278     int pos = conf_[output];
0279     if (pos >= (int)n_)
0280       std::copy(&deriv_.front() + n_ * (pos - n_), &deriv_.front() + n_ * (pos - n_ + 1), std::back_inserter(derivs));
0281     else {
0282       derivs.resize(n_);
0283       derivs[pos] = 1.;
0284     }
0285 
0286     return values_[pos];
0287   }
0288 
0289   template void MVAComputer::evalInternal(EvalContext &ctx) const;
0290   template void MVAComputer::evalInternal(DerivContext &ctx) const;
0291 
0292 }  // namespace PhysicsTools