HybridGenerator: Asyncronous + parallel event generation

This introduces:
* asyncronous event generation
* possibility for parallel event generation

This is useful for:
* hiding latency (IO) of certain generators
* decoupling the actual work from the call sequence into HybridGenerator
* collaboration from multiple clones of the same generator to generate
  a certain number of events

The implementation relies on tbb::task_arena and input/output queues
for decoupling the task_arena from the HybridGenerator thread.

Small adjustments to seeding of Pythia8 in order
to avoid same seeds in multiple parallel Pythia instances.

Author: sawenzel

Generators/include/Generators/Generator.h

@@ -141,6 +141,11 @@ class Generator : public FairGenerator
   /** lorentz boost data members **/
   Double_t mBoost;
 
+  // a unique generator instance counter
+  // this can be used to make sure no two generator instances have the same seed etc.
+  static std::atomic<int> InstanceCounter;
+  int mThisInstanceID = 0;
+
  private:
   void updateSubGeneratorInformation(o2::dataformats::MCEventHeader* header) const;
 

Generators/include/Generators/GeneratorHybrid.h

@@ -39,6 +39,12 @@
 #include <rapidjson/writer.h>
 #include "TBufferJSON.h"
 
+#include <tbb/concurrent_queue.h>
+#include <tbb/task_arena.h>
+#include <iostream>
+#include <thread>
+#include <atomic>
+
 namespace o2
 {
 namespace eventgen
@@ -84,6 +90,24 @@ class GeneratorHybrid : public Generator
   int mseqCounter = 0;
   int mCurrentFraction = 0;
   int mIndex = 0;
+
+  // Create a task arena with a specified number of threads
+  std::thread mTBBTaskPoolRunner;
+  tbb::concurrent_bounded_queue<int> mInputTaskQueue;
+  std::vector<tbb::concurrent_bounded_queue<int>> mResultQueue;
+  tbb::task_arena mTaskArena;
+  std::atomic<bool> mStopFlag;
+  bool mIsInitialized = false;
+
+  enum class GenMode {
+    kSeq,
+    kParallel
+  };
+
+  // hybrid gen operation mode - should be either 'sequential' or 'parallel'
+  // parallel means that we have clones of the same generator collaborating on event generation
+  // sequential means that events will be produced in the order given by fractions; async processing is still happening
+  GenMode mGenerationMode = GenMode::kSeq; //!
 };
 
 } // namespace eventgen

Generators/include/Generators/GeneratorHybridParam.h

@@ -31,6 +31,7 @@ namespace eventgen
 struct GeneratorHybridParam : public o2::conf::ConfigurableParamHelper<GeneratorHybridParam> {
   std::string configFile = ""; // JSON configuration file for the generators
   bool randomize = false;      // randomize the order of the generators, if not generator using fractions
+  int num_workers = 1;         // number of threads available for asyn/parallel event generation
   O2ParamDef(GeneratorHybridParam, "GeneratorHybrid");
 };
 

Generators/include/Generators/GeneratorPythia8.h

@@ -287,6 +287,9 @@ class GeneratorPythia8 : public Generator
                                // Value of -1 means unitialized; 0 will be time-dependent and values >1 <= MAX_SEED concrete reproducible seeding
   Pythia8GenConfig mGenConfig; // configuration object
 
+  static std::atomic<int> Pythia8InstanceCounter;
+  int mThisPythia8InstanceID = 0;
+
   constexpr static long MAX_SEED = 900000000;
 
   ClassDefOverride(GeneratorPythia8, 1);

Generators/src/Generator.cxx

@@ -28,13 +28,17 @@ namespace o2
 namespace eventgen
 {
 
+std::atomic<int> Generator::InstanceCounter{0};
+
 /*****************************************************************/
 /*****************************************************************/
 
 Generator::Generator() : FairGenerator("ALICEo2", "ALICEo2 Generator"),
                          mBoost(0.)
 {
   /** default constructor **/
+  mThisInstanceID = Generator::InstanceCounter;
+  Generator::InstanceCounter++;
 }
 
 /*****************************************************************/
@@ -43,6 +47,8 @@ Generator::Generator(const Char_t* name, const Char_t* title) : FairGenerator(na
                                                                 mBoost(0.)
 {
   /** constructor **/
+  mThisInstanceID = Generator::InstanceCounter;
+  Generator::InstanceCounter++;
 }
 
 /*****************************************************************/

Generators/src/GeneratorHybrid.cxx

@@ -13,6 +13,10 @@
 #include <fairlogger/Logger.h>
 #include <algorithm>
 
+#include <tbb/concurrent_queue.h>
+#include <tbb/task_arena.h>
+#include <tbb/parallel_for.h>
+
 namespace o2
 {
 namespace eventgen
@@ -34,7 +38,7 @@ GeneratorHybrid::GeneratorHybrid(const std::string& inputgens)
     }
   }
   int index = 0;
-  if (!mRandomize) {
+  if (!(mRandomize || mGenerationMode == GenMode::kParallel)) {
     if (mFractions.size() != mInputGens.size()) {
       LOG(fatal) << "Number of fractions does not match the number of generators";
       return;
@@ -159,54 +163,116 @@ Bool_t GeneratorHybrid::Init()
   } else {
     LOG(info) << "Generators will be used in sequence, following provided fractions";
   }
+
+  if (mGenerationMode == GenMode::kParallel) {
+    // in parallel mode we just use one queue --> collaboration
+    mResultQueue.resize(1);
+  } else {
+    // in sequential mode we have one queue per generator
+    mResultQueue.resize(gens.size());
+  }
+  // Create a task arena with a specified number of threads
+  mTaskArena.initialize(GeneratorHybridParam::Instance().num_workers);
+
+  // the process task function actually calls event generation
+  // when it is done it notifies the outside world by pushing it's index into an appropriate queue
+  // This should be a lambda, which can be given at TaskPool creation time
+  auto process_generator_task = [this](int task) {
+    LOG(info) << "Starting eventgen for task " << task;
+    auto& generator = gens[task];
+    generator->clearParticles();
+    generator->generateEvent();
+    generator->importParticles();
+    if (mGenerationMode == GenMode::kParallel) {
+      mResultQueue[0].push(task);
+    } else {
+      mResultQueue[task].push(task);
+    }
+  };
+
+  // fundamental tbb thread-worker function
+  auto worker_function = [this, process_generator_task]() {
+    while (!mStopFlag) {
+      int task;
+      if (mInputTaskQueue.try_pop(task)) {
+        process_generator_task(task); // Process the task
+      } else {
+        std::this_thread::sleep_for(std::chrono::milliseconds(10)); // Wait if no task
+      }
+    }
+    std::cout << "Worker thread exiting on thread " << std::this_thread::get_id() << std::endl;
+  };
+
+  // let the TBB task system run in it's own thread
+  mTBBTaskPoolRunner = std::thread([this, worker_function]() { mTaskArena.execute([&]() { tbb::parallel_for(0, mTaskArena.max_concurrency(), [&](int) { worker_function(); }); }); });
+
+  mTBBTaskPoolRunner.detach();
+
+  // let's initially push initial generation tasks for each event generator
+  for (size_t genindex = 0; genindex < gens.size(); ++genindex) {
+    mInputTaskQueue.push(genindex);
+  }
+  mIsInitialized = true;
   return Generator::Init();
 }
 
-Bool_t GeneratorHybrid::generateEvent()
+bool GeneratorHybrid::generateEvent()
 {
-  // Order randomisation or sequence of generators
-  // following provided fractions. If not available generators will be used sequentially
-  if (mRandomize) {
-    if (mRngFractions.size() != 0) {
-      // Generate number between 0 and 1
-      float rnum = gRandom->Rndm();
-      // Find generator index
-      for (int k = 0; k < mRngFractions.size(); k++) {
-        if (rnum <= mRngFractions[k]) {
-          mIndex = k;
-          break;
+  if (!mIsInitialized) {
+    Init();
+  }
+  if (mGenerationMode == GenMode::kParallel) {
+    mIndex = -1;           // this means any index is welcome
+    notifySubGenerator(0); // we shouldn't distinguish the sub-gen ids
+  } else {
+    // Order randomisation or sequence of generators
+    // following provided fractions, if not generators are used in proper sequence
+    // Order randomisation or sequence of generators
+    // following provided fractions. If not available generators will be used sequentially
+    if (mRandomize) {
+      if (mRngFractions.size() != 0) {
+        // Generate number between 0 and 1
+        float rnum = gRandom->Rndm();
+        // Find generator index
+        for (int k = 0; k < mRngFractions.size(); k++) {
+          if (rnum <= mRngFractions[k]) {
+            mIndex = k;
+            break;
+          }
         }
+      } else {
+        mIndex = gRandom->Integer(mGens.size());
       }
     } else {
-      mIndex = gRandom->Integer(mGens.size());
-    }
-  } else {
-    while (mFractions[mCurrentFraction] == 0 || mseqCounter == mFractions[mCurrentFraction]) {
-      if (mFractions[mCurrentFraction] != 0) {
-        mseqCounter = 0;
+      while (mFractions[mCurrentFraction] == 0 || mseqCounter == mFractions[mCurrentFraction]) {
+        if (mFractions[mCurrentFraction] != 0) {
+          mseqCounter = 0;
+        }
+        mCurrentFraction = (mCurrentFraction + 1) % mFractions.size();
       }
-      mCurrentFraction = (mCurrentFraction + 1) % mFractions.size();
+      mIndex = mCurrentFraction;
     }
-    mIndex = mCurrentFraction;
+    notifySubGenerator(mIndex);
   }
-  if (mConfigs[mIndex].compare("") == 0) {
-    LOG(info) << "GeneratorHybrid: generating event with generator " << mGens[mIndex];
-  } else {
-    LOG(info) << "GeneratorHybrid: generating event with generator " << mConfigs[mIndex];
-  }
-  gens[mIndex]->clearParticles(); // clear container of this class
-  gens[mIndex]->generateEvent();
-  // notify the sub event generator
-  notifySubGenerator(mIndex);
-  mseqCounter++;
   return true;
 }
 
-Bool_t GeneratorHybrid::importParticles()
+bool GeneratorHybrid::importParticles()
 {
-  mParticles.clear(); // clear container of mother class
-  gens[mIndex]->importParticles();
-  std::copy(gens[mIndex]->getParticles().begin(), gens[mIndex]->getParticles().end(), std::back_insert_iterator(mParticles));
+  int genIndex = -1;
+  if (mIndex == -1) {
+    // this means parallel mode ---> we have a common queue
+    mResultQueue[0].pop(genIndex);
+  } else {
+    // need to pop from a particular queue
+    mResultQueue[mIndex].pop(genIndex);
+  }
+
+  // at this moment the mIndex-th generator is ready to be used
+  std::copy(gens[genIndex]->getParticles().begin(), gens[genIndex]->getParticles().end(), std::back_insert_iterator(mParticles));
+
+  // we initiate the next async event generation stage for this generator
+  mInputTaskQueue.push(genIndex);
 
   // we need to fix particles statuses --> need to enforce this on the importParticles level of individual generators
   for (auto& p : mParticles) {
@@ -215,6 +281,7 @@ Bool_t GeneratorHybrid::importParticles()
     p.SetBit(ParticleStatus::kToBeDone, true);
   }
 
+  mseqCounter++;
   return true;
 }
 
@@ -244,6 +311,21 @@ Bool_t GeneratorHybrid::parseJSON(const std::string& path)
     return false;
   }
 
+  // check if there is a mode field
+  if (doc.HasMember("mode")) {
+    const auto& mode = doc["mode"].GetString();
+    if (mode == "sequential") {
+      // events are generated in the order given by fractions or random weight
+      mGenerationMode = GenMode::kSeq;
+    }
+    if (mode == std::string("parallel")) {
+      // events are generated fully in parallel and the order will be random
+      // this is mainly for event pool generation or mono-type generators
+      mGenerationMode = GenMode::kParallel;
+      LOG(info) << "Setting mode to parallel";
+    }
+  }
+
   // Put the generator names in mInputGens
   if (doc.HasMember("generators")) {
     const auto& gens = doc["generators"];

Generators/src/GeneratorPythia8.cxx

@@ -42,6 +42,8 @@ namespace o2
 namespace eventgen
 {
 
+std::atomic<int> GeneratorPythia8::Pythia8InstanceCounter;
+
 /*****************************************************************/
 /*****************************************************************/
 
@@ -57,6 +59,8 @@ GeneratorPythia8::GeneratorPythia8() : GeneratorPythia8(GeneratorPythia8Param::I
 GeneratorPythia8::GeneratorPythia8(Pythia8GenConfig const& config) : Generator("ALICEo2", "ALICEo2 Pythia8 Generator")
 {
   /** constructor **/
+  mThisPythia8InstanceID = GeneratorPythia8::Pythia8InstanceCounter;
+  GeneratorPythia8::Pythia8InstanceCounter++;
 
   mInterface = reinterpret_cast<void*>(&mPythia);
   mInterfaceName = "pythia8";
@@ -116,7 +120,15 @@ void GeneratorPythia8::seedGenerator()
     // Otherwise will seed the generator with the state of
     // TRandom::GetSeed. This is the seed that is influenced from
     // SimConfig --seed command line options options.
-    seed = (gRandom->TRandom::GetSeed() % (MAX_SEED + 1));
+    seed = gRandom->TRandom::GetSeed(); // this uses the "original" seed
+    // we advance the seed by one so that the next Pythia8 generator gets a different value
+    if (mThisPythia8InstanceID > 0) {
+      gRandom->Rndm();
+      LOG(info) << "Multiple Pythia8 generator instances detected .. automatically adjusting seed further to avoid overlap ";
+      seed = seed ^ gRandom->GetSeed(); // this uses the "current" seed
+    }
+    // apply max seed cuttof
+    seed = seed % (MAX_SEED + 1);
     LOG(info) << "GeneratorPythia8: Using random seed from gRandom % 900000001: " << seed;
   }
   mPythia.readString("Random:setSeed on");