Dependency injection with FermiBreakUp and Configurations

CMakeLists.txt

@@ -1,7 +1,7 @@
 cmake_minimum_required(VERSION 3.21)
 project(FermiBreakUp)
 
-set(CMAKE_CXX_STANDARD 11)
+set(CMAKE_CXX_STANDARD 17)
 
 find_package(CLHEP REQUIRED)
 
@@ -31,11 +31,13 @@ set(Fragments
         MyFermiBreakUp/Fragments/He5FermiFragment.cpp)
 
 set(Sources
+        MyFermiBreakUp/VFermiBreakUp.cpp
         MyFermiBreakUp/IntegerPartition.cpp
         MyFermiBreakUp/FermiParticle.cpp
         MyFermiBreakUp/FermiSplit.cpp
         MyFermiBreakUp/FermiFragmentPool.cpp
         MyFermiBreakUp/FermiConfigurations.cpp
+        MyFermiBreakUp/CachedFermiConfigurations.cpp
         MyFermiBreakUp/FermiBreakUp.cpp)
 
 set(Includes ${Sources} ${Fragments} ${Utilities})
@@ -60,7 +62,8 @@ set(Tests
         MyFermiBreakUp/Test/PartitionTest.cpp
         MyFermiBreakUp/Test/TableTests.cpp
         MyFermiBreakUp/Test/BreakUpTest.cpp
-        MyFermiBreakUp/Test/LRUTest.cpp)
+        MyFermiBreakUp/Test/LRUTest.cpp
+        MyFermiBreakUp/Test/CachedBreakUpTest.cpp)
 
 add_executable(CTest ${Tests} ${Includes})
 

MyFermiBreakUp/CachedFermiConfigurations.cpp

@@ -0,0 +1,64 @@
+//
+// Created by Artem Novikov on 24.05.2023.
+//
+
+#include "CachedFermiConfigurations.h"
+#include "ConfigurationProperties.h"
+#include "Randomizer.h"
+
+CachedFermiConfigurations::CachedFermiConfigurations(NucleiData nuclei_data, FermiFloat total_energy) {
+  CachedFermiConfigurations::GenerateSplits(nuclei_data, total_energy);
+}
+
+VFermiConfigurations& CachedFermiConfigurations::GenerateSplits(NucleiData nuclei_data, FermiFloat total_energy) {
+  auto max_fragments_count = FermiUInt(nuclei_data.mass_number);
+  if (nuclei_data != last_nuclei_) {
+    last_nuclei_ = nuclei_data;
+    cached_configurations_.clear();
+
+    for (uint32_t particle_count = 2; particle_count <= max_fragments_count; particle_count++) {
+      for (auto& split : FermiSplit(nuclei_data, particle_count)) {
+        cached_configurations_.emplace_back(std::move(split));  /// split is moved!
+      }
+    }
+  }
+
+  configurations_.clear();
+  weights_.clear();
+
+  FermiFloat total_weight = 0;
+  for(size_t i = 0; i < cached_configurations_.size(); ++i) {
+    auto split_weight = ConfigurationProperties::DecayProbability(cached_configurations_[i], nuclei_data.mass_number, total_energy);
+
+    if (split_weight != 0) {
+      total_weight += split_weight;
+
+      weights_.push_back(split_weight);
+      configurations_.emplace_back(i);
+    }
+  }
+
+  std::transform(weights_.begin(), weights_.end(),
+                 weights_.begin(), std::bind(std::divides<FermiFloat>(), std::placeholders::_1, total_weight));
+
+  return *this;
+}
+
+std::optional<FragmentVector> CachedFermiConfigurations::ChooseSplit() {
+  if (configurations_.empty()) {
+    return {};
+  }
+
+  FermiFloat wheel_result = Randomizer::UniformRealDistribution();
+  FermiFloat accumulated_weight = 0;
+
+  for (size_t i = 0; i < weights_.size(); ++i) {
+    accumulated_weight += weights_[i];
+
+    if (accumulated_weight >= wheel_result) {
+      return cached_configurations_[configurations_[i]];
+    }
+  }
+
+  throw std::runtime_error("No split chosen, something went wrong!");
+}

MyFermiBreakUp/CachedFermiConfigurations.h

@@ -0,0 +1,29 @@
+//
+// Created by Artem Novikov on 24.05.2023.
+//
+
+#ifndef FERMIBREAKUP_MYFERMIBREAKUP_CACHEDFERMICONFIGURATIONS_H_
+#define FERMIBREAKUP_MYFERMIBREAKUP_CACHEDFERMICONFIGURATIONS_H_
+
+#include "FermiSplit.h"
+#include "VFermiConfigurations.h"
+
+class CachedFermiConfigurations : public VFermiConfigurations {
+ public:
+  CachedFermiConfigurations() = default;
+
+  CachedFermiConfigurations(NucleiData nuclei_data, FermiFloat total_energy);
+
+  VFermiConfigurations& GenerateSplits(NucleiData nuclei_data, FermiFloat total_energy) override;
+
+  std::optional<FragmentVector> ChooseSplit() override;
+
+ private:
+  std::vector<size_t> configurations_;
+  std::vector<FermiFloat> weights_;
+
+  std::vector<FragmentVector> cached_configurations_;
+  NucleiData last_nuclei_{};
+};
+
+#endif //FERMIBREAKUP_MYFERMIBREAKUP_CACHEDFERMICONFIGURATIONS_H_

MyFermiBreakUp/FermiBreakUp.cpp

@@ -4,9 +4,13 @@
 
 #include "FermiBreakUp.h"
 #include "DataTypes.h"
-#include "PhaseDecay/FermiPhaseSpaceDecay.h"
 
-ParticleSplit FermiBreakUp::BreakItUp(const FermiParticle& nucleus) {
+FermiBreakUp::FermiBreakUp() : fermi_configurations_(DefaultConfigurations()) {}
+
+FermiBreakUp::FermiBreakUp(std::unique_ptr<VFermiConfigurations>&& configurations)
+    : fermi_configurations_(std::move(configurations)) {}
+
+ParticleVector FermiBreakUp::BreakItUp(const FermiParticle& nucleus) {
   /// CHECK that Excitation Energy > 0
   if (nucleus.GetExcitationEnergy() < 0) {
     return {nucleus};
@@ -16,37 +20,14 @@ ParticleSplit FermiBreakUp::BreakItUp(const FermiParticle& nucleus) {
   FermiFloat total_energy = nucleus.GetMomentum().m();
 
   /// Split the nucleus
-  FermiConfigurations configurations(nucleus.GetMassNumber(), nucleus.GetChargeNumber(), total_energy);
-  if (!configurations.IsSplitPossible()) {
+  auto fragment_split = fermi_configurations_->GenerateSplits(nucleus.GetNucleiData(), total_energy).ChooseSplit();
+  if (!fragment_split.has_value()) {
     return {nucleus};
   }
 
-  return ConvertToParticles(nucleus, configurations.ChooseSplit());
+  return ConvertToParticles(nucleus, fragment_split.value());
 }
 
-ParticleSplit FermiBreakUp::ConvertToParticles(const FermiParticle& source_nucleus, const FragmentSplit& split) {
-  ParticleSplit particle_split;
-  particle_split.reserve(2 * split.size());
-
-  std::vector<FermiFloat> split_masses;
-  split_masses.reserve(split.size());
-  for (auto fragment_ptr : split) {
-    split_masses.push_back(fragment_ptr->GetTotalEnergy());
-  }
-
-  FermiPhaseSpaceDecay phase_sampler;
-  std::vector<LorentzVector> particles_momentum = phase_sampler.CalculateDecay(source_nucleus.GetMomentum(), split_masses);
-
-  Vector3 boost_vector = source_nucleus.GetMomentum().boostVector();
-
-  /// Go back to the Lab Frame
-  for (size_t fragment_idx = 0; fragment_idx < split.size(); ++fragment_idx) {
-    ParticleVector fragment_particles = split[fragment_idx]->GetFragment(
-        particles_momentum[fragment_idx].boost(boost_vector));
-
-    particle_split.insert(particle_split.end(), std::make_move_iterator(fragment_particles.begin()),
-              std::make_move_iterator(fragment_particles.end()));
-  }
-
-  return particle_split;
+std::unique_ptr<VFermiConfigurations> FermiBreakUp::DefaultConfigurations() {
+  return std::make_unique<FermiConfigurations>();
 }

MyFermiBreakUp/FermiBreakUp.h

@@ -7,19 +7,23 @@
 
 #include "FermiParticle.h"
 #include "FermiConfigurations.h"
+#include "VFermiBreakUp.h"
+#include "VFermiConfigurations.h"
 
-using ParticleSplit = std::vector<FermiParticle>;
-
-class FermiBreakUp {
+class FermiBreakUp : public VFermiBreakUp {
  public:
-  FermiBreakUp() = delete;
+  FermiBreakUp();
+
+  FermiBreakUp(std::unique_ptr<VFermiConfigurations>&& configurations);
 
-  static ParticleSplit BreakItUp(const FermiParticle& nucleus);
+  ParticleVector BreakItUp(const FermiParticle& nucleus) override;
 
   ~FermiBreakUp() = default;
 
  private:
-  static ParticleSplit ConvertToParticles(const FermiParticle& source_nucleus, const FragmentSplit& split);
+  static std::unique_ptr<VFermiConfigurations> DefaultConfigurations();
+
+  std::unique_ptr<VFermiConfigurations> fermi_configurations_;
 };
 
 #endif //FERMIBREAKUP_MYFERMIBREAKUP_FERMIBREAKUP_H_

MyFermiBreakUp/FermiConfigurations.cpp

@@ -8,17 +8,23 @@
 #include "ConfigurationProperties.h"
 #include "Randomizer.h"
 
-FermiConfigurations::FermiConfigurations(MassNumber mass_number, ChargeNumber charge_number, FermiFloat total_energy) {
+FermiConfigurations::FermiConfigurations(NucleiData nuclei_data, FermiFloat total_energy) {
+  FermiConfigurations::GenerateSplits(nuclei_data, total_energy);
+}
+
+VFermiConfigurations& FermiConfigurations::GenerateSplits(NucleiData nuclei_data, FermiFloat total_energy) {
+  configurations_.clear();
+  weights_.clear();
+
   /// let's split nucleus into 2,...,A fragments
   FermiFloat total_weight = 0;
-  auto max_fragments_count = FermiUInt(mass_number);
+  auto max_fragments_count = FermiUInt(nuclei_data.mass_number);
 
   for (uint32_t particle_count = 2; particle_count <= max_fragments_count; particle_count++) {
     /// initialize configuration for k fragments
-    /// TODO cache splits
-    for (auto& split : FermiSplit(mass_number, charge_number, particle_count)) {
+    for (auto& split : FermiSplit(nuclei_data, particle_count)) {
       /// non-normalized statistical weight for given channel with k fragments
-      auto split_weight = ConfigurationProperties::DecayProbability(split, mass_number, total_energy);
+      auto split_weight = ConfigurationProperties::DecayProbability(split, nuclei_data.mass_number, total_energy);
       if (split_weight != 0) {
         total_weight += split_weight;
 
@@ -30,16 +36,14 @@ FermiConfigurations::FermiConfigurations(MassNumber mass_number, ChargeNumber ch
 
   /// let's normalize statistical weights of channels
   std::transform(weights_.begin(), weights_.end(),
-                 weights_.begin(), std::bind(std::divides<FermiFloat>(), std::placeholders::_1, total_weight));
-}
+                 weights_.begin(), std::bind(std::divides(), std::placeholders::_1, total_weight));
 
-bool FermiConfigurations::IsSplitPossible() const {
-  return !configurations_.empty();
+  return *this;
 }
 
-const FragmentSplit& FermiConfigurations::ChooseSplit() const {
-  if (!IsSplitPossible()) {
-    throw std::runtime_error("No possible splits!");
+std::optional<FragmentVector> FermiConfigurations::ChooseSplit() {
+  if (configurations_.empty()) {
+    return {};
   }
 
   FermiFloat wheel_result = Randomizer::UniformRealDistribution();

MyFermiBreakUp/FermiConfigurations.h

@@ -8,18 +8,20 @@
 #include <vector>
 
 #include "FermiSplit.h"
+#include "VFermiConfigurations.h"
 
-class FermiConfigurations {
+class FermiConfigurations : public VFermiConfigurations {
  public:
-  FermiConfigurations(MassNumber mass_number, ChargeNumber charge_number, FermiFloat total_energy);
+  FermiConfigurations() = default;
 
-  bool IsSplitPossible() const;
+  FermiConfigurations(NucleiData nuclei_data, FermiFloat total_energy);
 
-  const FragmentSplit& ChooseSplit() const;
+  VFermiConfigurations& GenerateSplits(NucleiData nuclei_data, FermiFloat total_energy) override;
 
- private:
+  std::optional<FragmentVector> ChooseSplit() override;
 
-  std::vector<FragmentSplit> configurations_;
+ private:
+  std::vector<FragmentVector> configurations_;
   std::vector<FermiFloat> weights_;
 };
 

MyFermiBreakUp/FermiFragmentPool.cpp

@@ -135,7 +135,7 @@ size_t FermiFragmentPool::Count(MassNumber mass_number, ChargeNumber charge_numb
   return Count({mass_number, charge_number});
 }
 
-size_t FermiFragmentPool::Count(const NucleiData& nuclei) const {
+size_t FermiFragmentPool::Count(NucleiData nuclei) const {
   return fragments_pool_->count(nuclei);
 }
 
@@ -148,6 +148,6 @@ FermiFragmentPool::RangeIterators FermiFragmentPool::GetFragments(MassNumber mas
   return GetFragments({mass_number, charge_number});
 }
 
-FermiFragmentPool::RangeIterators FermiFragmentPool::GetFragments(const NucleiData& nuclei) const {
+FermiFragmentPool::RangeIterators FermiFragmentPool::GetFragments(NucleiData nuclei) const {
   return fragments_pool_->equal_range(nuclei);
 }

MyFermiBreakUp/FermiFragmentPool.h

@@ -21,11 +21,11 @@ class FermiFragmentPool {
 
   size_t Count(MassNumber mass_number, ChargeNumber charge_number) const;
 
-  size_t Count(const NucleiData& nuclei) const;
+  size_t Count(NucleiData nuclei) const;
 
   RangeIterators GetFragments(MassNumber mass_number, ChargeNumber charge_number) const;
 
-  RangeIterators GetFragments(const NucleiData& nuclei) const;
+  RangeIterators GetFragments(NucleiData nuclei) const;
 
  private:
   void AddFragment(const FermiFragment& fragment) const;

MyFermiBreakUp/FermiParticle.cpp

@@ -14,6 +14,11 @@ FermiParticle::FermiParticle(MassNumber mass_number, ChargeNumber charge_number,
   CalculateExcitationEnergy();
 }
 
+NucleiData FermiParticle::GetNucleiData() const {
+  return {mass_number_, charge_number_};
+}
+
+
 MassNumber FermiParticle::GetMassNumber() const {
   return mass_number_;
 }

MyFermiBreakUp/FermiParticle.h

@@ -23,6 +23,8 @@ class FermiParticle {
 
   FermiParticle(MassNumber mass_number, ChargeNumber charge_number, const LorentzVector& momentum);
 
+  NucleiData GetNucleiData() const;
+
   MassNumber GetMassNumber() const;
 
   ChargeNumber GetChargeNumber() const;