File indexing completed on 2024-04-06 12:23:38
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017 #include <cstdlib>
0018 #include <algorithm>
0019 #include <iterator>
0020 #include <vector>
0021 #include <cmath>
0022
0023 #include "FWCore/Utilities/interface/Exception.h"
0024
0025 #include "PhysicsTools/MVAComputer/interface/VarProcessor.h"
0026 #include "PhysicsTools/MVAComputer/interface/Calibration.h"
0027
0028 using namespace PhysicsTools;
0029
0030 namespace {
0031
0032 class ProcMLP : public VarProcessor {
0033 public:
0034 typedef VarProcessor::Registry::Registry<ProcMLP, Calibration::ProcMLP> Registry;
0035
0036 ProcMLP(const char *name, const Calibration::ProcMLP *calib, const MVAComputer *computer);
0037 ~ProcMLP() override {}
0038
0039 void configure(ConfIterator iter, unsigned int n) override;
0040 void eval(ValueIterator iter, unsigned int n) const override;
0041 std::vector<double> deriv(ValueIterator iter, unsigned int n) const override;
0042
0043 private:
0044 struct Layer {
0045 Layer(const Calibration::ProcMLP::Layer &calib);
0046 Layer(const Layer &orig)
0047 : inputs(orig.inputs), neurons(orig.neurons), coeffs(orig.coeffs), sigmoid(orig.sigmoid) {}
0048
0049 unsigned int inputs;
0050 unsigned int neurons;
0051 std::vector<double> coeffs;
0052 bool sigmoid;
0053 };
0054
0055 std::vector<Layer> layers;
0056 unsigned int maxTmp;
0057 };
0058
0059 ProcMLP::Registry registry("ProcMLP");
0060
0061 ProcMLP::Layer::Layer(const Calibration::ProcMLP::Layer &calib)
0062 : inputs(calib.first.front().second.size()), neurons(calib.first.size()), sigmoid(calib.second) {
0063 typedef Calibration::ProcMLP::Neuron Neuron;
0064
0065 coeffs.resize(neurons * (inputs + 1));
0066 std::vector<double>::iterator inserter = coeffs.begin();
0067
0068 for (std::vector<Neuron>::const_iterator iter = calib.first.begin(); iter != calib.first.end(); iter++) {
0069 *inserter++ = iter->first;
0070
0071 if (iter->second.size() != inputs)
0072 throw cms::Exception("ProcMLPInput") << "ProcMLP neuron layer inconsistent." << std::endl;
0073
0074 inserter = std::copy(iter->second.begin(), iter->second.end(), inserter);
0075 }
0076 }
0077
0078 ProcMLP::ProcMLP(const char *name, const Calibration::ProcMLP *calib, const MVAComputer *computer)
0079 : VarProcessor(name, calib, computer), maxTmp(0) {
0080 std::copy(calib->layers.begin(), calib->layers.end(), std::back_inserter(layers));
0081
0082 for (unsigned int i = 0; i < layers.size(); i++) {
0083 maxTmp = std::max<unsigned int>(maxTmp, layers[i].neurons);
0084 if (i > 0 && layers[i - 1].neurons != layers[i].inputs)
0085 throw cms::Exception("ProcMLPInput") << "ProcMLP neuron layers do not connect "
0086 "properly."
0087 << std::endl;
0088 }
0089 }
0090
0091 void ProcMLP::configure(ConfIterator iter, unsigned int n) {
0092 if (n != layers.front().inputs)
0093 return;
0094
0095 for (unsigned int i = 0; i < n; i++)
0096 iter++(Variable::FLAG_NONE);
0097
0098 for (unsigned int i = 0; i < layers.back().neurons; i++)
0099 iter << Variable::FLAG_NONE;
0100 }
0101
0102 void ProcMLP::eval(ValueIterator iter, unsigned int n) const {
0103 double *tmp = (double *)alloca(2 * maxTmp * sizeof(double));
0104 bool flip = false;
0105
0106 for (double *pos = tmp; iter; iter++, pos++)
0107 *pos = *iter;
0108
0109 double *output = nullptr;
0110 for (std::vector<Layer>::const_iterator layer = layers.begin(); layer != layers.end(); layer++, flip = !flip) {
0111 const double *input = &tmp[flip ? maxTmp : 0];
0112 output = &tmp[flip ? 0 : maxTmp];
0113 std::vector<double>::const_iterator coeff = layer->coeffs.begin();
0114 for (unsigned int i = 0; i < layer->neurons; i++) {
0115 double sum = *coeff++;
0116 for (unsigned int j = 0; j < layer->inputs; j++)
0117 sum += input[j] * *coeff++;
0118 if (layer->sigmoid)
0119 sum = 1.0 / (std::exp(-sum) + 1.0);
0120 *output++ = sum;
0121 }
0122 }
0123
0124 for (const double *pos = &tmp[flip ? maxTmp : 0]; pos < output; pos++)
0125 iter(*pos);
0126 }
0127
0128 std::vector<double> ProcMLP::deriv(ValueIterator iter, unsigned int n) const {
0129 std::vector<double> prevValues, nextValues;
0130 std::vector<double> prevMatrix, nextMatrix;
0131
0132 while (iter)
0133 nextValues.push_back(*iter++);
0134
0135 unsigned int size = nextValues.size();
0136 nextMatrix.resize(size * size);
0137 for (unsigned int i = 0; i < size; i++)
0138 nextMatrix[i * size + i] = 1.;
0139
0140 for (std::vector<Layer>::const_iterator layer = layers.begin(); layer != layers.end(); layer++) {
0141 prevValues.clear();
0142 std::swap(prevValues, nextValues);
0143 prevMatrix.clear();
0144 std::swap(prevMatrix, nextMatrix);
0145
0146 std::vector<double>::const_iterator coeff = layer->coeffs.begin();
0147 for (unsigned int i = 0; i < layer->neurons; i++) {
0148 double sum = *coeff++;
0149 for (unsigned int j = 0; j < layer->inputs; j++)
0150 sum += prevValues[j] * *coeff++;
0151
0152 double deriv;
0153 if (layer->sigmoid) {
0154 double e = std::exp(-sum);
0155 sum = 1.0 / (e + 1.0);
0156 deriv = 1.0 / (e + 1.0 / e + 2.0);
0157 } else
0158 deriv = 1.0;
0159
0160 nextValues.push_back(sum);
0161
0162 for (unsigned int k = 0; k < size; k++) {
0163 sum = 0.0;
0164 coeff -= layer->inputs;
0165 for (unsigned int j = 0; j < layer->inputs; j++)
0166 sum += prevMatrix[j * size + k] * *coeff++;
0167 nextMatrix.push_back(sum * deriv);
0168 }
0169 }
0170 }
0171
0172 return nextMatrix;
0173 }
0174
0175 }