174 {
175
176 Weightfile weightfile;
177 std::string custom_weightfile = weightfile.generateFileName();
178 std::string custom_steeringfile = weightfile.generateFileName();
179
180 uint64_t numberOfFeatures = training_data.getNumberOfFeatures();
181 uint64_t numberOfSpectators = training_data.getNumberOfSpectators();
182 uint64_t numberOfEvents = training_data.getNumberOfEvents();
183
184 if (m_specific_options.m_training_fraction <= 0.0 or m_specific_options.m_training_fraction > 1.0) {
185 B2ERROR("Please provide a positive training fraction");
186 throw std::runtime_error("Please provide a training fraction between (0.0,1.0]");
187 }
188
189 auto numberOfTrainingEvents = static_cast<uint64_t>(numberOfEvents * 100 * m_specific_options.m_training_fraction);
190 numberOfTrainingEvents = numberOfTrainingEvents / 100 + (numberOfTrainingEvents % 100 != 0);
191 auto numberOfValidationEvents = numberOfEvents - numberOfTrainingEvents;
192
193 uint64_t batch_size = m_specific_options.m_mini_batch_size;
194 if (batch_size == 0) {
195 batch_size = numberOfTrainingEvents;
196 }
197
198 if (batch_size > numberOfTrainingEvents) {
199 B2WARNING("Mini batch size (" << batch_size << ") is larger than the number of training events (" << numberOfTrainingEvents << ")"\
200 " The batch size has been set equal to the number of training events.");
201 batch_size = numberOfTrainingEvents;
202 };
203
204 auto X = std::unique_ptr<float[]>(new float[batch_size * numberOfFeatures]);
205 auto S = std::unique_ptr<float[]>(new float[batch_size * numberOfSpectators]);
206 auto y = std::unique_ptr<float[]>(new float[batch_size]);
207 auto w = std::unique_ptr<float[]>(new float[batch_size]);
208 npy_intp dimensions_X[2] = {static_cast<npy_intp>(batch_size), static_cast<npy_intp>(numberOfFeatures)};
209 npy_intp dimensions_S[2] = {static_cast<npy_intp>(batch_size), static_cast<npy_intp>(numberOfSpectators)};
210 npy_intp dimensions_y[2] = {static_cast<npy_intp>(batch_size), 1};
211 npy_intp dimensions_w[2] = {static_cast<npy_intp>(batch_size), 1};
212
213 auto X_v = std::unique_ptr<float[]>(new float[numberOfValidationEvents * numberOfFeatures]);
214 auto S_v = std::unique_ptr<float[]>(new float[numberOfValidationEvents * numberOfSpectators]);
215 auto y_v = std::unique_ptr<float[]>(new float[numberOfValidationEvents]);
216 auto w_v = std::unique_ptr<float[]>(new float[numberOfValidationEvents]);
217 npy_intp dimensions_X_v[2] = {static_cast<npy_intp>(numberOfValidationEvents), static_cast<npy_intp>(numberOfFeatures)};
218 npy_intp dimensions_S_v[2] = {static_cast<npy_intp>(numberOfValidationEvents), static_cast<npy_intp>(numberOfSpectators)};
219 npy_intp dimensions_y_v[2] = {static_cast<npy_intp>(numberOfValidationEvents), 1};
220 npy_intp dimensions_w_v[2] = {static_cast<npy_intp>(numberOfValidationEvents), 1};
221
222 std::string steering_file_source_code;
223 if (m_specific_options.m_steering_file != "") {
224 steering_file_source_code = loadPythonFileAsString(m_specific_options.m_steering_file);
225 }
226
227 std::vector<float> means(numberOfFeatures, 0.0);
228 std::vector<float> stds(numberOfFeatures, 0.0);
229
230 if (m_specific_options.m_normalize) {
231
232
233 auto weights = training_data.getWeights();
234 for (uint64_t iFeature = 0; iFeature < numberOfFeatures; ++iFeature) {
235 double wSum = 0.0;
236 double mean = 0.0;
237 double running_std = 0.0;
238 auto feature = training_data.getFeature(iFeature);
239 for (uint64_t i = 0; i < weights.size(); ++i) {
240 wSum += weights[i];
241 double meanOld = mean;
242 mean += (weights[i] / wSum) * (feature[i] - meanOld);
243 running_std += weights[i] * (feature[i] - meanOld) * (feature[i] - mean);
244 }
245 means[iFeature] = mean;
246 stds[iFeature] = std::sqrt(running_std / (wSum - 1));
247 }
248 }
249
250 try {
251
252 auto json = boost::python::import("json");
253 auto builtins = boost::python::import("builtins");
254 auto inspect = boost::python::import("inspect");
255
256
257
258 boost::python::object type = boost::python::import("types");
259
260
261 boost::uuids::random_generator uuid_gen;
262 std::string unique_mva_module_name = "unique_module_name" + boost::uuids::to_string(uuid_gen());
263 boost::python::object unique_mva_module = type.attr("ModuleType")(unique_mva_module_name.c_str());
264
265
266 auto framework = boost::python::import((std::string("basf2_mva_python_interface.") + m_specific_options.m_framework).c_str());
267 auto framework_file = framework.attr("__file__");
268 boost::python::extract<std::string> extractor(framework_file);
269 std::string framework_filename = extractor();
270 auto framework_file_source_code = loadPythonFileAsString(framework_filename);
271 builtins.attr("exec")(framework_file_source_code.c_str(), boost::python::object(unique_mva_module.attr("__dict__")));
272
273 builtins.attr("exec")(steering_file_source_code.c_str(), boost::python::object(unique_mva_module.attr("__dict__")));
274
275
276 auto parameters = json.attr("loads")(m_specific_options.m_config.c_str());
277 auto model = unique_mva_module.attr("get_model")(numberOfFeatures, numberOfSpectators,
278 numberOfEvents, m_specific_options.m_training_fraction, parameters);
279
280
281 for (uint64_t iEvent = 0; iEvent < numberOfValidationEvents; ++iEvent) {
282 training_data.loadEvent(iEvent);
283 if (m_specific_options.m_normalize) {
284 for (uint64_t iFeature = 0; iFeature < numberOfFeatures; ++iFeature)
285 X_v[iEvent * numberOfFeatures + iFeature] = (training_data.m_input[iFeature] - means[iFeature]) / stds[iFeature];
286 } else {
287 for (uint64_t iFeature = 0; iFeature < numberOfFeatures; ++iFeature)
288 X_v[iEvent * numberOfFeatures + iFeature] = training_data.m_input[iFeature];
289 }
290 for (uint64_t iSpectator = 0; iSpectator < numberOfSpectators; ++iSpectator)
291 S_v[iEvent * numberOfSpectators + iSpectator] = training_data.m_spectators[iSpectator];
292 y_v[iEvent] = training_data.m_target;
293 w_v[iEvent] = training_data.m_weight;
294 }
295
296 auto ndarray_X_v = boost::python::handle<>(PyArray_SimpleNewFromData(2, dimensions_X_v, NPY_FLOAT32, X_v.get()));
297 auto ndarray_S_v = boost::python::handle<>(PyArray_SimpleNewFromData(2, dimensions_S_v, NPY_FLOAT32, S_v.get()));
298 auto ndarray_y_v = boost::python::handle<>(PyArray_SimpleNewFromData(2, dimensions_y_v, NPY_FLOAT32, y_v.get()));
299 auto ndarray_w_v = boost::python::handle<>(PyArray_SimpleNewFromData(2, dimensions_w_v, NPY_FLOAT32, w_v.get()));
300
301 uint64_t nBatches = std::floor(numberOfTrainingEvents / batch_size);
302
303 auto state = unique_mva_module.attr("begin_fit")(model, ndarray_X_v, ndarray_S_v, ndarray_y_v, ndarray_w_v, nBatches);
304
305 bool continue_loop = true;
306
307 std::vector<uint64_t> iteration_index_vector(numberOfTrainingEvents);
308 std::iota(std::begin(iteration_index_vector), std::end(iteration_index_vector), 0);
309
310 for (uint64_t iIteration = 0; (iIteration < m_specific_options.m_nIterations or m_specific_options.m_nIterations == 0)
311 and continue_loop; ++iIteration) {
312
313
314 if (iIteration > 0) std::shuffle(std::begin(iteration_index_vector), std::end(iteration_index_vector), TRandomWrapper());
315
316 for (uint64_t iBatch = 0; iBatch < nBatches and continue_loop; ++iBatch) {
317
318
319
320 Py_BEGIN_ALLOW_THREADS;
321 for (uint64_t iEvent = 0; iEvent < batch_size; ++iEvent) {
322 training_data.loadEvent(iteration_index_vector.at(iEvent + iBatch * batch_size) + numberOfValidationEvents);
323 if (m_specific_options.m_normalize) {
324 for (uint64_t iFeature = 0; iFeature < numberOfFeatures; ++iFeature)
325 X[iEvent * numberOfFeatures + iFeature] = (training_data.m_input[iFeature] - means[iFeature]) / stds[iFeature];
326 } else {
327 for (uint64_t iFeature = 0; iFeature < numberOfFeatures; ++iFeature)
328 X[iEvent * numberOfFeatures + iFeature] = training_data.m_input[iFeature];
329 }
330 for (uint64_t iSpectator = 0; iSpectator < numberOfSpectators; ++iSpectator)
331 S[iEvent * numberOfSpectators + iSpectator] = training_data.m_spectators[iSpectator];
332 y[iEvent] = training_data.m_target;
333 w[iEvent] = training_data.m_weight;
334 }
335
336
337 Py_END_ALLOW_THREADS;
338
339
340 auto ndarray_X = boost::python::handle<>(PyArray_SimpleNewFromData(2, dimensions_X, NPY_FLOAT32, X.get()));
341 auto ndarray_S = boost::python::handle<>(PyArray_SimpleNewFromData(2, dimensions_S, NPY_FLOAT32, S.get()));
342 auto ndarray_y = boost::python::handle<>(PyArray_SimpleNewFromData(2, dimensions_y, NPY_FLOAT32, y.get()));
343 auto ndarray_w = boost::python::handle<>(PyArray_SimpleNewFromData(2, dimensions_w, NPY_FLOAT32, w.get()));
344
345 auto r = unique_mva_module.attr("partial_fit")(state, ndarray_X, ndarray_S, ndarray_y,
346 ndarray_w, iIteration, iBatch);
347 boost::python::extract<bool> proxy(r);
348 if (proxy.check())
349 continue_loop = static_cast<bool>(proxy);
350 }
351 }
352
353 auto result = unique_mva_module.attr("end_fit")(state);
354
355 auto pickle = boost::python::import("pickle");
356 auto file = builtins.attr("open")(custom_weightfile.c_str(), "wb");
357 pickle.attr("dump")(result, file);
358
359 auto steeringfile = builtins.attr("open")(custom_steeringfile.c_str(), "wb");
360 pickle.attr("dump")(steering_file_source_code.c_str(), steeringfile);
361
362 auto importances = unique_mva_module.attr("feature_importance")(state);
363 if (len(importances) == 0) {
364 B2INFO("Python method returned empty feature importance. There won't be any information about the feature importance in the weightfile.");
365 } else if (numberOfFeatures != static_cast<uint64_t>(len(importances))) {
366 B2WARNING("Python method didn't return the correct number of importance value. I ignore the importances");
367 } else {
368 std::map<std::string, float> feature_importances;
369 for (uint64_t iFeature = 0; iFeature < numberOfFeatures; ++iFeature) {
370 boost::python::extract<float> proxy(importances[iFeature]);
371 if (proxy.check()) {
372 feature_importances[m_general_options.m_variables[iFeature]] = static_cast<float>(proxy);
373 } else {
374 B2WARNING("Failed to convert importance output of the method to a float, using 0 instead");
375 feature_importances[m_general_options.m_variables[iFeature]] = 0.0;
376 }
377 }
378 weightfile.addFeatureImportance(feature_importances);
379 }
380
381 } catch (...) {
382 PyErr_Print();
383 PyErr_Clear();
384 B2ERROR("Failed calling train in PythonTeacher");
385 throw std::runtime_error(std::string("Failed calling train in PythonTeacher"));
386 }
387
388 weightfile.addOptions(m_general_options);
389 weightfile.addOptions(m_specific_options);
390 weightfile.addFile("Python_Weightfile", custom_weightfile);
391 weightfile.addFile("Python_Steeringfile", custom_steeringfile);
392 weightfile.addSignalFraction(training_data.getSignalFraction());
393 if (m_specific_options.m_normalize) {
394 weightfile.addVector("Python_Means", means);
395 weightfile.addVector("Python_Stds", stds);
396 }
397
398 return weightfile;
399
400 }