Merge pull request #285 from OpenGATE/update_gaga_garf

update garf + gaga (WIP do not merge)

.github/workflows/main.yml

@@ -17,10 +17,10 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        os: [ubuntu-latest, macos-11, windows-latest]
+        os: [ubuntu-latest, macos-latest, windows-latest]
         python-version: [3.8, 3.9, '3.10', '3.11']
         exclude:
-          - os: macos-11
+          - os: macos-latest
             python-version: '3.11'
 
     env:
@@ -36,7 +36,7 @@ jobs:
         if [[ ${{ matrix.os }} == "ubuntu-latest" ]]; then
           export GIT_SSL_NO_VERIFY=1
           git submodule update --init --recursive
-        elif [[ ${{ matrix.os }} == "macos-11" ]]; then
+        elif [[ ${{ matrix.os }} == "macos-latest" ]]; then
           export GIT_SSL_NO_VERIFY=1
           git submodule update --init --recursive
         else
@@ -56,7 +56,7 @@ jobs:
             varOS=`cat /etc/os-release | grep "VERSION=" | grep -oP '(?<=\").*?(?=\")'`
             varOS=($varOS)
             echo "release=${varOS[0]}" >> $GITHUB_OUTPUT
-        elif [[ ${{ matrix.os }} == "macos-11" ]]; then
+        elif [[ ${{ matrix.os }} == "macos-latest" ]]; then
             varOS=`sw_vers | grep "ProductVersion:"`
             varOS="${varOS#*:}"
             echo "release=${varOS:1}" >> $GITHUB_OUTPUT
@@ -106,22 +106,22 @@ jobs:
               mv dist_opengate/* dist/
           fi
     - uses: conda-incubator/setup-miniconda@v2
-      if: (matrix.os == 'macos-11') || (matrix.os == 'windows-latest')
+      if: (matrix.os == 'macos-latest') || (matrix.os == 'windows-latest')
       with:
           auto-update-conda: true
           activate-environment: opengate_core
           python-version: ${{ matrix.python-version }}
     - name: Set up Homebrew
-      if: matrix.os == 'macos-11'
+      if: matrix.os == 'macos-latest'
       id: set-up-homebrew
       uses: Homebrew/actions/setup-homebrew@master
     - name: Create opengate_core Wheel Mac
-      if: matrix.os == 'macos-11'
+      if: matrix.os == 'macos-latest'
       shell: bash -l {0}
       run: |
           brew update
-          rm -rf /usr/local/bin/python3.1*-config /usr/local/bin/2to3-3.1* /usr/local/bin/idle3.1* /usr/local/bin/pydoc3.1* /usr/local/bin/python3.1*
-          rm -rf /usr/local/bin/python3-config /usr/local/bin/2to3 /usr/local/bin/idle3 /usr/local/bin/pydoc3 /usr/local/bin/python3
+          #rm -rf /usr/local/bin/python3.1*-config /usr/local/bin/2to3-3.1* /usr/local/bin/idle3.1* /usr/local/bin/pydoc3.1* /usr/local/bin/python3.1*
+          #rm -rf /usr/local/bin/python3-config /usr/local/bin/2to3 /usr/local/bin/idle3 /usr/local/bin/pydoc3 /usr/local/bin/python3
           brew install --force --verbose --overwrite \
                        ccache \
                        fftw \
@@ -363,8 +363,7 @@ jobs:
             pip install torch
           fi
           pip install SimpleITK
-          pip install gaga_phsp
-          pip install garf
+          pip install gaga_phsp==0.7.1
           pip install dist/opengate_core-*-${PYTHONFOLDER}-${OSNAME}*_${PLATFORM}64.whl
           pip install dist/opengate-*.whl
           export GIT_SSL_NO_VERIFY=1

core/opengate_core/opengate_lib/GateARFActor.cpp

@@ -16,6 +16,7 @@ GateARFActor::GateARFActor(py::dict &user_info) : GateVActor(user_info, true) {
   fActions.insert("EndOfRunAction");
   // User option: batch size
   fBatchSize = DictGetInt(user_info, "batch_size");
+  fKeepNegativeSide = DictGetBool(user_info, "flip_plane");
 }
 
 void GateARFActor::SetARFFunction(ARFFunctionType &f) { fApply = f; }
@@ -36,17 +37,31 @@ void GateARFActor::EndOfRunAction(const G4Run * /*run*/) {
     l.fPositionY.clear();
     l.fDirectionX.clear();
     l.fDirectionY.clear();
+    l.fDirectionZ.clear();
     l.fCurrentNumberOfHits = 0;
   }
 }
 
 void GateARFActor::SteppingAction(G4Step *step) {
   auto &l = fThreadLocalData.Get();
 
+  // get direction and transform to local
+  auto *pre = step->GetPreStepPoint();
+  auto dir = pre->GetMomentumDirection();
+  dir = pre->GetTouchable()->GetHistory()->GetTopTransform().TransformAxis(dir);
+
+  // which side of the plane ?
+  if (!fKeepNegativeSide && dir[2] < 0)
+    return;
+  if (fKeepNegativeSide && dir[2] > 0)
+    return;
+
   l.fCurrentNumberOfHits++;
+  l.fDirectionX.push_back(dir[0]);
+  l.fDirectionY.push_back(dir[1]);
+  // l.fDirectionZ.push_back(dir[2]); // not used
 
   // get energy
-  auto *pre = step->GetPreStepPoint();
   l.fEnergy.push_back(pre->GetKineticEnergy());
 
   // get position and transform to local
@@ -56,12 +71,6 @@ void GateARFActor::SteppingAction(G4Step *step) {
   l.fPositionX.push_back(pos[0]);
   l.fPositionY.push_back(pos[1]);
 
-  // get direction and transform to local
-  auto dir = pre->GetMomentumDirection();
-  dir = pre->GetTouchable()->GetHistory()->GetTopTransform().TransformAxis(dir);
-  l.fDirectionX.push_back(dir[0]);
-  l.fDirectionY.push_back(dir[1]);
-
   // trigger the "apply" (ARF) if the number of hits in the batch is reached
   if (l.fCurrentNumberOfHits >= fBatchSize) {
     fApply(this);
@@ -70,6 +79,7 @@ void GateARFActor::SteppingAction(G4Step *step) {
     l.fPositionY.clear();
     l.fDirectionX.clear();
     l.fDirectionY.clear();
+    l.fDirectionZ.clear();
     l.fCurrentNumberOfHits = 0;
   }
 }
@@ -101,3 +111,7 @@ std::vector<double> GateARFActor::GetDirectionX() const {
 std::vector<double> GateARFActor::GetDirectionY() const {
   return fThreadLocalData.Get().fDirectionY;
 }
+
+std::vector<double> GateARFActor::GetDirectionZ() const {
+  return fThreadLocalData.Get().fDirectionZ;
+}

core/opengate_core/opengate_lib/GateARFActor.h

@@ -43,6 +43,8 @@ class GateARFActor : public GateVActor {
 
   std::vector<double> GetDirectionY() const;
 
+  std::vector<double> GetDirectionZ() const;
+
   // Main function called every step in attached volume
   void SteppingAction(G4Step *) override;
 
@@ -52,6 +54,7 @@ class GateARFActor : public GateVActor {
 protected:
   int fBatchSize;
   ARFFunctionType fApply;
+  bool fKeepNegativeSide;
 
   // For MT, all threads local variables are gathered here
   struct threadLocalT {
@@ -60,6 +63,7 @@ class GateARFActor : public GateVActor {
     std::vector<double> fPositionY;
     std::vector<double> fDirectionX;
     std::vector<double> fDirectionY;
+    std::vector<double> fDirectionZ;
     // number of particle hitting the detector
     int fCurrentNumberOfHits;
     // Current run id (to detect if run has changed)

core/opengate_core/opengate_lib/GateGANSource.cpp

@@ -238,6 +238,7 @@ G4ThreeVector GateGANSource::GeneratePrimariesPosition() {
     position =
         G4ThreeVector(fPositionX[fCurrentIndex], fPositionY[fCurrentIndex],
                       fPositionZ[fCurrentIndex]);
+    position = fLocalRotation * position + fLocalTranslation; // FIXME
     // move position according to mother volume
     auto &l = fThreadLocalData.Get();
     position = l.fGlobalRotation * position + l.fGlobalTranslation;

core/opengate_core/opengate_lib/GateHelpersGeometry.cpp

@@ -24,7 +24,11 @@ void ComputeTransformationFromVolumeToWorld(const std::string &phys_volume_name,
     auto rot = phys->GetObjectRotationValue();
     rotation = rot * rotation;
     translation = rot * translation + tr;
-    name = phys->GetMotherLogical()->GetName();
+    // Warning, the world can be a parallel world
+    if (phys->GetMotherLogical() == nullptr)
+      name = "world";
+    else
+      name = phys->GetMotherLogical()->GetName();
   }
 }
 

core/opengate_core/opengate_lib/pyGateARFActor.cpp

@@ -24,5 +24,6 @@ void init_GateARFActor(py::module &m) {
       .def("GetPositionX", &GateARFActor::GetPositionX)
       .def("GetPositionY", &GateARFActor::GetPositionY)
       .def("GetDirectionX", &GateARFActor::GetDirectionX)
-      .def("GetDirectionY", &GateARFActor::GetDirectionY);
+      .def("GetDirectionY", &GateARFActor::GetDirectionY)
+      .def("GetDirectionZ", &GateARFActor::GetDirectionZ);
 }

opengate/actors/arfactors.py

@@ -30,7 +30,7 @@ def import_garf():
     from packaging import version
 
     garf_version = pkg_resources.get_distribution("garf").version
-    garf_minimal_version = "2.4"
+    garf_minimal_version = "2.5"
     if version.parse(garf_version) < version.parse(garf_minimal_version):
         fatal(
             "The minimal version of garf is not correct. You should install at least the version "
@@ -105,13 +105,16 @@ def set_default_user_info(user_info):
         user_info.verbose_batch = False
         user_info.output = ""
         user_info.enable_hit_slice = False
+        user_info.flip_plane = False
         # Can be cpu / auto / gpu
         user_info.gpu_mode = "auto"
 
     def __init__(self, user_info):
         ActorBase.__init__(self, user_info)
         g4.GateARFActor.__init__(self, user_info.__dict__)
         # import module
+        self.debug_nb_hits_before = None
+        self.debug_nb_hits = 0
         self.garf = import_garf()
         if self.garf is None:
             print("Cannot run GANSource")
@@ -152,13 +155,17 @@ def initialize(self, volume_engine=None):
         self.ActorInitialize()
         self.SetARFFunction(self.apply)
         self.user_info.output_image = None
+        self.debug_nb_hits_before = 0
+        self.debug_nb_hits = 0
 
         # load the pth file
         self.nn, self.model = self.garf.load_nn(self.pth_filename, verbose=False)
         p = self.param
         p.batch_size = int(float(self.user_info.batch_size))
 
-        # size and spacing (2D)
+        # size and spacing (2D) (force to float)
+        self.user_info.image_spacing[0] = float(self.user_info.image_spacing[0])
+        self.user_info.image_spacing[1] = float(self.user_info.image_spacing[1])
         p.image_size = self.user_info.image_size
         p.image_spacing = self.user_info.image_spacing
         p.distance_to_crystal = self.user_info.distance_to_crystal
@@ -223,13 +230,15 @@ def apply_with_lock(self, actor):
 
         # build the data
         x = np.column_stack((px, py, theta, phi, energy))
+        self.debug_nb_hits_before += len(x)
 
         # apply the neural network
         if self.user_info.verbose_batch:
             print(
                 f"Apply ARF to {energy.shape[0]} hits (device = {self.model_data['current_gpu_mode']})"
             )
-        ax = x[:, 2:5]  # two angles and energy
+
+        ax = x[:, 2:5]  # two angles and energy # FIXME index ?
         w = self.garf.nn_predict(self.model, self.nn["model_data"], ax)
 
         # positions
@@ -242,20 +251,17 @@ def apply_with_lock(self, actor):
         coord = np.around(coord).astype(int)
         v = coord[:, 0]
         u = coord[:, 1]
-        u, v, w_pred = self.garf.remove_out_of_image_boundaries(u, v, w, p.image_size)
+        u, v, w_pred = self.garf.remove_out_of_image_boundaries2(
+            u, v, w, self.user_info.image_size
+        )
 
         # do nothing if there is no hit in the image
         if u.shape[0] != 0:
-            temp = np.zeros(p.image_size, dtype=np.float64)
-            temp = self.garf.image_from_coordinates(temp, u, v, w_pred)
-            # add to previous, at the correct slice location
-            # the slice is : current_ene_window + run_id * nb_ene_windows
             run_id = actor.GetCurrentRunId()
-            # self.simulation_engine_wr().g4_RunManager.GetCurrentRun().GetRunID()
             s = p.nb_ene * run_id
-            self.output_image[s : s + p.nb_ene] = (
-                self.output_image[s : s + p.nb_ene] + temp
-            )
+            img = self.output_image[s : s + p.nb_ene]
+            self.garf.image_from_coordinates_add(img, u, v, w_pred)
+            self.debug_nb_hits += u.shape[0]
 
     def EndSimulationAction(self):
         g4.GateARFActor.EndSimulationAction(self)
@@ -265,7 +271,7 @@ def EndSimulationAction(self):
         # Should we keep the first slice (with all hits) ?
         if not self.user_info.enable_hit_slice:
             self.output_image = self.output_image[1:, :, :]
-            self.param.image_size[1] = self.param.image_size[1] - 1
+            self.param.image_size[0] = self.param.image_size[0] - 1
 
         # convert to itk image
         self.output_image = itk.image_from_array(self.output_image)
@@ -294,3 +300,7 @@ def EndSimulationAction(self):
             itk.imwrite(
                 self.output_image, ensure_filename_is_str(self.user_info.output)
             )
+
+        # debug
+        # print(f"{self.debug_nb_hits_before=}")
+        # print(f"{self.debug_nb_hits=}")

opengate/actors/doseactors.py

@@ -508,7 +508,7 @@ def EndSimulationAction(self):
             )
             itk.imwrite(self.py_LETd_image, ensure_filename_is_str(fPath))
 
-            # for parrallel computation we need to provide both outputs
+            # for parallel computation we need to provide both outputs
             if self.user_info.separate_output:
                 fPath = fPath.replace(".mhd", "_numerator.mhd")
                 itk.imwrite(self.py_numerator_image, ensure_filename_is_str(fPath))

opengate/actors/miscactors.py

@@ -4,9 +4,8 @@
 import uproot
 import numpy as np
 import time
-
+import platform
 import opengate_core as g4
-
 from .base import ActorBase
 from ..exception import fatal
 from ..geometry.utility import rot_np_as_g4, vec_np_as_g4, vec_g4_as_np
@@ -49,6 +48,7 @@ def __init__(self, user_info=None):
         self.counts.stop_time = 0
         self.counts.init = 0
         self.counts.track_types = {}
+        self.nb_thread = 1
 
     @property
     def pps(self):
@@ -71,14 +71,6 @@ def sps(self):
             return self.counts.step_count / self.counts.duration * sec
         return 0
 
-    @property
-    def nb_thread(self):
-        if self.simulation is not None:
-            thread = self.simulation.number_of_threads
-        else:
-            thread = "?"
-        return thread
-
     @property
     def simu_start_time(self):
         if not self.simulation is None:
@@ -109,16 +101,22 @@ def __str__(self):
             f"PPS       {self.pps:.0f}\n"
             f"TPS       {self.tps:.0f}\n"
             f"SPS       {self.sps:.0f}\n"
-            f"start     {self.counts.start_time}\n"
-            f"stop      {self.counts.stop_time}\n"
+            f"Start     {self.counts.start_time}\n"
+            f"Stop      {self.counts.stop_time}\n"
             f'Sim start {g4.G4BestUnit(self.simu_start_time, "Time")}\n'
             f'Sim end   {g4.G4BestUnit(self.simu_end_time, "Time")}\n'
-            f"Threads   {self.nb_thread}"
+            f"Threads   {self.nb_thread}\n"
+            f"Arch      {platform.system()}\n"
+            f"Python    {platform.python_version()}"
         )
         if self.user_info.track_types_flag:
             s += f"\n" f"Track types: {self.counts.track_types}"
         return s
 
+    def StartSimulationAction(self):
+        g4.GateSimulationStatisticsActor.StartSimulationAction(self)
+        self.nb_thread = self.simulation.user_info.number_of_threads
+
     def EndSimulationAction(self):
         g4.GateSimulationStatisticsActor.EndSimulationAction(self)
         self.counts = Box(self.GetCounts())
@@ -162,6 +160,8 @@ def write(self, filename):
         s += f"# SPS (Step per sec)         = {self.sps:.0f}\n"
         s += f"# Threads                    = {self.nb_thread}\n"
         s += f"# Date                       = {datetime.now()}\n"
+        s += f"# Arch                       = {platform.system()}\n"
+        s += f"# Python                     = {platform.python_version()}\n"
         if self.user_info.track_types_flag:
             s += f"# Track types:\n"
             for t in self.counts.track_types:

opengate/bin/dose_rate

@@ -35,10 +35,10 @@ def go(json_param, output_folder):
     sim = create_simulation(param)
 
     # run
-    output = sim.start()
+    sim.run()
 
     # print results at the end
-    stats = output.get_actor("Stats")
+    stats = sim.output.get_actor("Stats")
     print(stats)
     stats.write(param.output_folder / "stats.txt")
     print(f"Output in {param.output_folder}")