Test dose actor variants mt

core/opengate_core/opengate_core.cpp

@@ -249,6 +249,8 @@ void init_GateKineticEnergyFilter(py::module &);
 
 void init_GateDoseActor(py::module &m);
 
+void init_GateDoseSpeedTestActor(py::module &m);
+
 void init_GateFluenceActor(py::module &m);
 
 void init_GateLETActor(py::module &m);
@@ -477,6 +479,7 @@ PYBIND11_MODULE(opengate_core, m) {
   init_GateEventAction(m);
   init_GateTrackingAction(m);
   init_GateDoseActor(m);
+  init_GateDoseSpeedTestActor(m);
   init_GateFluenceActor(m);
   init_GateLETActor(m);
   init_GateSimulationStatisticsActor(m);

core/opengate_core/opengate_lib/GateDoseSpeedTestActor.cpp

@@ -0,0 +1,199 @@
+/* --------------------------------------------------
+   Copyright (C): OpenGATE Collaboration
+   This software is distributed under the terms
+   of the GNU Lesser General  Public Licence (LGPL)
+   See LICENSE.md for further details
+   ------------------------------------ -------------- */
+
+#include "GateDoseSpeedTestActor.h"
+#include "G4Navigator.hh"
+#include "G4RandomTools.hh"
+#include "G4RunManager.hh"
+#include "GateHelpers.h"
+#include "GateHelpersDict.h"
+#include "GateHelpersImage.h"
+
+#include "G4Deuteron.hh"
+#include "G4Electron.hh"
+#include "G4EmCalculator.hh"
+#include "G4Gamma.hh"
+#include "G4MaterialTable.hh"
+#include "G4NistManager.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4ParticleTable.hh"
+#include "G4Positron.hh"
+#include "G4Proton.hh"
+#include "G4Threading.hh"
+#include "G4Types.hh"
+
+// Mutex that will be used by thread to write in the edep/dose image
+// TODO
+G4Mutex EndOfRunMutex = G4MUTEX_INITIALIZER;
+
+GateDoseSpeedTestActor::GateDoseSpeedTestActor(py::dict &user_info)
+    : GateVActor(user_info, true) {
+  // Create the image pointer
+  // (the size and allocation will be performed on the py side)
+  cpp_reference_image = ImageType::New();
+  // Action for this actor: during stepping
+  fActions.insert("SteppingAction");
+  fActions.insert("BeginOfRunAction");
+  fActions.insert("EndOfRunAction");
+  fActions.insert("EndSimulationAction");
+  fInitialTranslation = DictGetG4ThreeVector(user_info, "translation");
+  fstorageMethod = DictGetStr(user_info, "storage_method");
+  fcountWriteAttempts = DictGetBool(user_info, "count_write_attempts");
+}
+
+GateDoseSpeedTestActor::~GateDoseSpeedTestActor() {
+  delete deposit_vector_atomic_pointer;
+}
+
+void GateDoseSpeedTestActor::ActorInitialize() {
+  fNumberOfThreads = G4Threading::GetNumberOfRunningWorkerThreads();
+  if (fNumberOfThreads == 0) {
+    fNumberOfThreads = 1;
+  }
+}
+
+void GateDoseSpeedTestActor::PrepareStorage() {
+  ImageType::RegionType region =
+      cpp_reference_image->GetLargestPossibleRegion();
+  fimageSize = region.GetSize();
+  fnumberOfVoxels = fimageSize[0] * fimageSize[1] * fimageSize[2];
+
+  cpp_image = ImageType::New();
+  cpp_image->CopyInformation(cpp_reference_image);
+  RegionType region_cpp_image;
+  region_cpp_image.SetIndex(
+      cpp_reference_image->GetLargestPossibleRegion().GetIndex());
+  region_cpp_image.SetSize(
+      cpp_reference_image->GetLargestPossibleRegion().GetSize());
+  cpp_image->SetRegions(region_cpp_image);
+  cpp_image->Allocate();
+
+  if (fstorageMethod == "atomic") {
+    for (int i = 0; i < fnumberOfVoxels; i++) {
+      deposit_vector.emplace_back();
+    }
+  } else if (fstorageMethod == "standard") {
+    for (int i = 0; i < fnumberOfVoxels; i++) {
+      deposit_vector_standard.emplace_back(0);
+    }
+  } else if (fstorageMethod == "atomic_vec_pointer") {
+    delete deposit_vector_atomic_pointer;
+    deposit_vector_atomic_pointer =
+        new std::vector<std::atomic<double>>(fnumberOfVoxels);
+  }
+}
+
+void GateDoseSpeedTestActor::PrepareStorageLocal() {
+  auto &l = fThreadLocalData.Get();
+  l.deposit_vector_local.resize(fnumberOfVoxels);
+  std::fill(l.deposit_vector_local.begin(), l.deposit_vector_local.end(), 0.0);
+}
+
+void GateDoseSpeedTestActor::BeginOfRunAction(const G4Run *) {
+
+  if (fstorageMethod == "local") {
+    PrepareStorageLocal();
+  }
+
+  // Important ! The volume may have moved, so we re-attach each run
+  AttachImageToVolume<ImageType>(cpp_reference_image, fPhysicalVolumeName,
+                                 fInitialTranslation);
+  // compute volume of a dose voxel
+  auto sp = cpp_reference_image->GetSpacing();
+  fVoxelVolume = sp[0] * sp[1] * sp[2];
+}
+
+void GateDoseSpeedTestActor::SteppingAction(G4Step *step) {
+  auto position = step->GetPostStepPoint()->GetPosition();
+  auto touchable = step->GetPreStepPoint()->GetTouchable();
+  auto localPosition =
+      touchable->GetHistory()->GetTransform(0).TransformPoint(position);
+
+  // convert G4ThreeVector to itk PointType
+  ImageType::PointType point;
+  point[0] = localPosition[0];
+  point[1] = localPosition[1];
+  point[2] = localPosition[2];
+
+  // get pixel index
+  ImageType::IndexType index;
+  bool isInside =
+      cpp_reference_image->TransformPhysicalPointToIndex(point, index);
+
+  // set value
+  if (isInside) {
+    auto w = step->GetTrack()->GetWeight();
+    auto edep = step->GetTotalEnergyDeposit() / CLHEP::MeV * w;
+    int index_1d = (int)cpp_reference_image->ComputeOffset(index);
+    ftotalDepositWrites++;
+    if (fstorageMethod == "atomic") {
+      if (fcountWriteAttempts == true) {
+        ftotalReattemptsAtomicAdd += atomic_add_double_return_reattempts(
+            deposit_vector[index_1d].dep, edep);
+      } else {
+        atomic_add_double(deposit_vector[index_1d].dep, edep);
+      }
+    } else if (fstorageMethod == "standard") {
+      deposit_vector_standard[index_1d] += edep;
+    } else if (fstorageMethod == "local") {
+      auto &l = fThreadLocalData.Get();
+      l.deposit_vector_local[index_1d] += edep;
+    } else if (fstorageMethod == "atomic_vec_pointer") {
+      if (fcountWriteAttempts == true) {
+        ftotalReattemptsAtomicAdd += atomic_add_double_return_reattempts(
+            (*deposit_vector_atomic_pointer)[index_1d], edep);
+      } else {
+        atomic_add_double((*deposit_vector_atomic_pointer)[index_1d], edep);
+      }
+    }
+  } // else : outside the image
+}
+
+void GateDoseSpeedTestActor::EndSimulationAction() {
+  if (fstorageMethod == "atomic" || fstorageMethod == "standard" ||
+      fstorageMethod == "atomic_vec_pointer") {
+    PrepareOutput();
+  }
+}
+
+void GateDoseSpeedTestActor::EndOfRunAction(const G4Run *run) {
+  if (fstorageMethod == "local") {
+    WriteOutputToImageLocal();
+  }
+}
+
+void GateDoseSpeedTestActor::WriteOutputToImageLocal() {
+  G4AutoLock mutex(&EndOfRunMutex);
+
+  int index_1d;
+  auto &l = fThreadLocalData.Get();
+  ImageIteratorType it(cpp_image, cpp_image->GetLargestPossibleRegion());
+  it.GoToBegin();
+  while (!it.IsAtEnd()) {
+    index_1d = (int)cpp_reference_image->ComputeOffset(it.GetIndex());
+    ImageAddValue<ImageType>(cpp_image, it.GetIndex(),
+                             l.deposit_vector_local[index_1d]);
+    ++it;
+  }
+}
+
+void GateDoseSpeedTestActor::PrepareOutput() {
+  int index_1d;
+  ImageIteratorType it(cpp_image, cpp_image->GetLargestPossibleRegion());
+  it.GoToBegin();
+  while (!it.IsAtEnd()) {
+    index_1d = (int)cpp_reference_image->ComputeOffset(it.GetIndex());
+    if (fstorageMethod == "atomic") {
+      it.Set(deposit_vector[index_1d].dep);
+    } else if (fstorageMethod == "standard") {
+      it.Set(deposit_vector_standard[index_1d]);
+    } else if (fstorageMethod == "atomic_vec_pointer") {
+      it.Set((*deposit_vector_atomic_pointer)[index_1d]);
+    }
+    ++it;
+  }
+}

core/opengate_core/opengate_lib/GateDoseSpeedTestActor.h

@@ -0,0 +1,171 @@
+/* --------------------------------------------------
+   Copyright (C): OpenGATE Collaboration
+   This software is distributed under the terms
+   of the GNU Lesser General  Public Licence (LGPL)
+   See LICENSE.md for further details
+   -------------------------------------------------- */
+
+#ifndef GateDoseSpeedTestActor_h
+#define GateDoseSpeedTestActor_h
+
+#include "G4Cache.hh"
+#include "G4NistManager.hh"
+#include "G4VPrimitiveScorer.hh"
+#include "GateHelpers.h"
+#include "GateVActor.h"
+#include "itkImage.h"
+#include "itkImageRegionIterator.h"
+#include <atomic>
+#include <deque>
+#include <memory>
+#include <pybind11/stl.h>
+#include <stdlib.h>
+#include <vector>
+
+#define CACHELINESIZE 64
+
+namespace py = pybind11;
+
+struct dose_bin {
+  alignas(CACHELINESIZE) std::atomic<double> dep;
+  dose_bin() { dep = 0; }
+
+  dose_bin &operator+=(const double a) {
+    atomic_add_double(dep, a);
+    return *this;
+  }
+};
+
+// From https://github.com/zhourrr/aligned-memory-allocator
+// A minimal implementation of an allocator for C++ Standard Library, which
+// allocates aligned memory (specified by the alignment argument).
+template <typename T, std::size_t alignment> class AlignedAllocator {
+public:
+  using value_type = T;
+
+public:
+  // According to Microsoft's documentation, default ctor is not required
+  // by C++ Standard Library.
+  AlignedAllocator() noexcept {};
+
+  template <typename U>
+  AlignedAllocator(const AlignedAllocator<U, alignment> &other) noexcept {};
+
+  template <typename U>
+  inline bool
+  operator==(const AlignedAllocator<U, alignment> &other) const noexcept {
+    return true;
+  }
+
+  template <typename U>
+  inline bool
+  operator!=(const AlignedAllocator<U, alignment> &other) const noexcept {
+    return false;
+  }
+
+  template <typename U> struct rebind {
+    using other = AlignedAllocator<U, alignment>;
+  };
+
+  // STL containers call this function to allocate unintialized memory block to
+  // store (no more than n) elements of type T (value_type).
+  inline value_type *allocate(const std::size_t n) const {
+    auto size = n;
+    /*
+      If you wish, for some strange reason, that the size of allocated buffer is
+      also aligned to alignment, uncomment the following statement.
+
+      Note: this increases the size of underlying memory, but STL containers
+      still treat it as a memory block of size n, i.e., STL containers will not
+      put more than n elements into the returned memory.
+    */
+    // size = (n + alignment - 1) / alignment * alignment;
+    value_type *ptr;
+    auto ret = posix_memalign((void **)&ptr, alignment, sizeof(T) * size);
+    if (ret != 0)
+      throw std::bad_alloc();
+    return ptr;
+  };
+
+  // STL containers call this function to free a memory block beginning at a
+  // specified position.
+  inline void deallocate(value_type *const ptr, std::size_t n) const noexcept {
+    free(ptr);
+  }
+};
+
+class GateDoseSpeedTestActor : public GateVActor {
+
+public:
+  // Constructor
+  GateDoseSpeedTestActor(py::dict &user_info);
+
+  ~GateDoseSpeedTestActor();
+
+  // Main function called every step in attached volume
+  virtual void SteppingAction(G4Step *) override;
+
+  // Called every time a Run starts (all threads)
+  virtual void BeginOfRunAction(const G4Run *run) override;
+
+  virtual void EndOfRunAction(const G4Run *run) override;
+
+  virtual void EndSimulationAction() override;
+
+  // Called from the python-side when engines are initialized
+  virtual void ActorInitialize() override;
+
+  // pre-fill the vectors into which dose is written
+  virtual void PrepareStorage();
+
+  virtual void PrepareStorageLocal();
+
+  // prepare the dose images
+  void PrepareOutput();
+
+  void WriteOutputToImageLocal();
+
+  // Image type is 3D float by default
+  // TODO double precision required
+  using ImageType = itk::Image<double, 3>;
+  using ImageIteratorType = itk::ImageRegionIterator<ImageType>;
+  using RegionType = ImageType::RegionType;
+
+  // The image is accessible on py side (shared by all threads)
+  ImageType::Pointer cpp_reference_image;
+  ImageType::Pointer cpp_image;
+
+  using deposit_map_type = std::deque<dose_bin>;
+  //  using deposit_map_type = std::deque<dose_bin, AlignedAllocator<dose_bin,
+  //  CACHELINESIZE>>;
+  deposit_map_type deposit_vector;
+
+  std::vector<double> deposit_vector_standard;
+
+  std::vector<std::atomic<double>> *deposit_vector_atomic_pointer{};
+
+  std::string fPhysicalVolumeName;
+
+  int GetTotalReattemptsAtomicAdd() { return (int)ftotalReattemptsAtomicAdd; }
+
+  int GetTotalDepositWrites() { return (int)ftotalDepositWrites; }
+
+protected:
+  struct threadLocalT {
+    std::vector<double> deposit_vector_local;
+  };
+  G4Cache<threadLocalT> fThreadLocalData;
+
+private:
+  G4ThreeVector fInitialTranslation;
+  std::string fstorageMethod;
+  bool fcountWriteAttempts;
+  int fnumberOfVoxels;
+  double fVoxelVolume;
+  ImageType::SizeType fimageSize;
+  int fNumberOfThreads;
+  std::atomic<int> ftotalReattemptsAtomicAdd = 0;
+  std::atomic<int> ftotalDepositWrites = 0;
+};
+
+#endif // GateDoseSpeedTestActor_h