search

OpenGATE / opengate

Gate 10 (beta)
http://www.opengatecollaboration.org
GNU Lesser General Public License v3.0
43 stars 38 forks source link

Process dies simulating high activity over more than 10 s - std::out_of_range - The queue is empty #325

Open franky180 opened 8 months ago

franky180 commented 8 months ago

Simulating small amounts of activity (<30 MBq) over 100 s works fine with my current simulation setup. Increasing the activity to 60 MBq and 100 s results in an error message namely terminate called after throwing an instance of 'std::out_of_range' and leads to the multiprocessing manager Exception Exception: The queue is empty. The spawned process probably died.. Is there anything I can do do change this?

Error message:

See full error message 
(opengate_env) (base) user@test:~/opengate_env/lib/python3.11/site-packages/opengate/point source/src$ python GE_Healthcare_NM_670_Discovery_simulation_point_source_v2.3.py 
Importing opengate (thread 816115) ... done


New simulation with output path:  /home//opengate_env/lib/python3.11/site-packages/opengate/point source/src/../output/GE_NM_670_Tc99m_lehr_A=60 MBq_t=100 s_d=20 cm_bins=512_FWHM=13.3 keV_threads=8
Importing opengate (thread 816149) ... done
Importing opengate (thread 816186) ... done
Simulation: create MTRunManager with 8 threads
Simulation: initialize Geometry
Simulation: initialize Physics
Simulation: initialize Source
Simulation: initialize Actors
Simulation: initialize G4RunManager
Simulation: (no volumes overlap checking)
Simulation: apply G4 command '/tracking/verbose 0'
--------------------------------------------------------------------------------
Simulation: START (in a new process)
G4WT6 > 
-------- WWWW ------- G4Exception-START -------- WWWW -------
G4Exception : GeomNav1002
      issued by : G4Navigator::ComputeStep()
Stuck Track: potential geometry or navigation problem.
  Track stuck, not moving for 10 steps.
  Current  phys volume: 'spect_collimator_core'
   - at position : (160.9619480602623,55.90606160357735,-240.2)
     in direction: (0.5396197832743437,0.1897880022041721,-0.8202383822513425)
    (local position: (160.9619480602623,55.90606160357735,17.49999999999997))
    (local direction: (0.5396197832743437,0.1897880022041721,-0.8202383822513425)).
  Previous phys volume: 'spect_collimator_hole_param'

  Likely geometry overlap - else navigation problem !
   *** Trying to get *unstuck* using a push - expanding step to 1e-07 (mm) ...       Potential overlap in geometry !

*** This is just a warning message. ***
-------- WWWW -------- G4Exception-END --------- WWWW -------

terminate called after throwing an instance of 'std::out_of_range'
  what():  basic_string::substr: __pos (which is 1) > this->size() (which is 0)
(in /home//opengate_env/lib/python3.11/site-packages/opengate/processing.py line 30)
The queue is empty. The spawned process probably died.
Traceback (most recent call last):
  File "/home//opengate_env/lib/python3.11/site-packages/opengate/processing.py", line 28, in dispatch_to_subprocess
    return q.get(block=False)
           ^^^^^^^^^^^^^^^^^^
  File "", line 2, in get
  File "/usr/lib/python3.11/multiprocessing/managers.py", line 837, in _callmethod
    raise convert_to_error(kind, result)
_queue.Empty

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home//opengate_env/lib/python3.11/site-packages/opengate/point source/src/GE_Healthcare_NM_670_Discovery_simulation_point_source_v2.3.py", line 268, in 
    create_simulation(g_isotope, g_collimator_type, g_activity, g_runtime, g_distance, g_channels, g_num_of_threads, g_visu, g_FWHM, g_plot)
  File "/home//opengate_env/lib/python3.11/site-packages/opengate/point source/src/GE_Healthcare_NM_670_Discovery_simulation_point_source_v2.3.py", line 245, in create_simulation
    output = sim.run(start_new_process=True)
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home//opengate_env/lib/python3.11/site-packages/opengate/managers.py", line 1284, in run
    self.output = dispatch_to_subprocess(self._run_simulation_engine, True)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home//opengate_env/lib/python3.11/site-packages/opengate/processing.py", line 30, in dispatch_to_subprocess
    fatal("The queue is empty. The spawned process probably died.")
  File "/home//opengate_env/lib/python3.11/site-packages/opengate/exception.py", line 29, in fatal
    raise Exception(s)
Exception: The queue is empty. The spawned process probably died.

Code

See code 

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import opengate as gate
import opengate.contrib.spect.genm670 as gate_spect
import uproot
import pathlib
import matplotlib.pyplot as plt
import numpy as np

"""
 This version is set up to reproduce the error with high activities.
"""

def format_magnitude(value, unit="Bq"):
    """
    Formats a given numerical value with a specified unit using SI unit prefixes.

    Parameters:
    - value (float or int): The numerical value to be formatted.
    - unit (str, optional): The unit to be appended to the formatted value (default is "Bq").

    Returns:
    - str: A formatted string representing the value with the appropriate SI unit prefix, e.g., "1.2 kBq".

    Example:
    >>> format_magnitude(1200, "Bq")
    '1.2 kBq'
    >>> format_magnitude(4500000, "Bq")
    '4.5 MBq'
    >>> format_magnitude(0, "Bq")
    '0 Bq'
    >>> format_magnitude(2500, "eV")
    '2.5 keV'

    Notes:
    - SI unit prefixes (k, M, G, T, P, E, Z, Y) are used to scale the value for readability.
    - The function automatically determines the appropriate prefix based on the magnitude of the input value.
    - If the input value is zero, the function returns '0' followed by the specified unit.

    Reference:
    - SI unit prefixes: https://en.wikipedia.org/wiki/Metric_prefix
    """

    prefixes = ['', 'k', 'M', 'G', 'T', 'P', 'E', 'Z', 'Y']

    if value == 0:
        return f"0 {unit}"

    magnitude = int(value).bit_length() // 10

    if magnitude < 0:
        magnitude = 0
    elif magnitude > len(prefixes) - 1:
        magnitude = len(prefixes) - 1

    scaled_value = value / (1e3 ** magnitude)

    # Adjust formatting based on the last digit
    if scaled_value % 1 == 0:
        formatted_value = f"{scaled_value:.0f}"
    else:
        formatted_value = f"{scaled_value:.1f}"

    return f"{formatted_value} {prefixes[magnitude]}{unit}"

def create_simulation(g_isotope = "Tc99m", g_collimator_type = "lehr", g_activity = 10e6, g_runtime = 60, g_distance = 10, g_channels = 512, g_num_of_threads = 1, g_visu = False, g_FWHM = 9.5 * 140.511e3/100, g_plot = False):

    # ----------------------------------------
    # Setup of the Geant4 units.
    # ----------------------------------------
    m = gate.g4_units.m
    mm = gate.g4_units.mm
    keV = gate.g4_units.keV
    mm = gate.g4_units.mm
    Bq = gate.g4_units.Bq
    cm = gate.g4_units.cm
    eV = gate.g4_units.eV
    sec = gate.g4_units.second

    # ----------------------------------------
    # Setup the output path.
    # ----------------------------------------
    output_path = str(pathlib.Path(__file__).parent.resolve() / ".." / "output" / f"GE_NM_670_{g_isotope}_{g_collimator_type}_A={format_magnitude(g_activity, 'Bq')}_t={format_magnitude(g_runtime, 's')}_d={format_magnitude(g_distance, 'cm')}_bins={g_channels}_FWHM={format_magnitude(g_FWHM, 'eV')}_threads={g_num_of_threads}")
    print()
    print("***************")
    print("New simulation with output path: ", output_path)

    # ----------------------------------------
    # Setup the instance of the Simulation object sim.
    # ----------------------------------------
    sim = gate.Simulation()

    # ----------------------------------------
    # Set the main options.
    # ----------------------------------------
    sim.g4_verbose = False
    sim.visu = g_visu
    sim.visu_type = "vrml"
    sim.visu_filename = "geant4VisuFile.wrl"
    sim.number_of_threads = g_num_of_threads
    sim.check_volumes_overlap = False
    sim.random_engine = "MersenneTwister"
    sim.random_seed = "auto"

    sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option4"
    sim.physics_manager.enable_decay = False
    sim.physics_manager.set_production_cut(volume_name="world", particle_name="all", value=10 * mm)

    key_of_interest = "TotalEnergyDeposit" # Key for kinetic energy
    scalings = 1000  # Scaled to keV. For different energies has to be scaled differently.

    # ----------------------------------------
    # Load source data.
    # ----------------------------------------
    w, e = gate.sources.generic.get_rad_gamma_energy_spectrum(g_isotope)

    # ----------------------------------------
    # Setup of the world volume
    # ----------------------------------------
    world = sim.world
    world.size = [1.5 * m, 1.5 * m, 1.5 * m]
    sim.world.material = "G4_AIR"

    # ----------------------------------------
    # Add the point source of 1 cm diameter.
    # ----------------------------------------
    point_source = sim.add_source("GenericSource", f"{g_isotope}-Source")
    point_source.particle = "gamma"
    point_source.energy.type = "spectrum_lines"
    point_source.energy.spectrum_energy = e
    point_source.energy.spectrum_weight = w
    point_source.activity = g_activity * Bq / sim.number_of_threads
    point_source.position.type = "sphere"
    point_source.position.radius = 0.5 * cm
    point_source.position.translation = [0, 0, 0 * cm]
    point_source.direction.type = "iso"

    # ----------------------------------------
    # Add the GE Healthcare NM 670 Discovery spect head into the simulation.
    # ----------------------------------------
    spect, crystal = gate_spect.add_ge_nm67_spect_head(sim, "spect", collimator_type=g_collimator_type, debug=False)
    trn = spect.size[2] / 2 # Translation by half the head size (z-direction)
    spect.translation = [0, 0, -(g_distance * cm + trn)] # Translation by the distance between source and detector

    # ----------------------------------------
    # Define local spect production cuts
    # ----------------------------------------
    sim.physics_manager.set_production_cut(volume_name="spect", particle_name="gamma", value=0.1 * mm)
    sim.physics_manager.set_production_cut(volume_name="spect", particle_name="electron", value=0.1 * mm)
    sim.physics_manager.set_production_cut(volume_name="spect", particle_name="positron", value=0.1 * mm)

    # ----------------------------------------
    # Actors and Digitizers
    # ----------------------------------------

    # Hits collection digitizer
    hc = sim.add_actor("DigitizerHitsCollectionActor", "Hits")
    hc.output = output_path + str(".root")
    hc.mother = crystal.name
    hc.attributes = ["PostPosition", "TotalEnergyDeposit", "TrackVolumeCopyNo", "GlobalTime", "KineticEnergy", "ProcessDefinedStep", "PreStepUniqueVolumeID"]

    # Singles collection
    sc = sim.add_actor("DigitizerAdderActor", "Singles")
    sc.mother = crystal.name
    sc.input_digi_collection = "Hits"
    sc.policy = "EnergyWinnerPosition"
    sc.skip_attributes = ["KineticEnergy", "ProcessDefinedStep"]
    sc.output = hc.output

    # Blurr actor
    bc = sim.add_actor("DigitizerBlurringActor", "Singles_with_blur")
    bc.output = hc.output 
    bc.input_digi_collection = "Singles"
    bc.blur_attribute = "TotalEnergyDeposit"
    bc.blur_method = "Gaussian"
    bc.blur_fwhm = g_FWHM * eV

    # Add efficiancy actor
    ea = sim.add_actor("DigitizerEfficiencyActor", "Efficiency")
    ea.mother = hc.mother 
    ea.input_digi_collection = "Singles_with_blur"
    ea.output = hc.output
    ea.efficiency = 0.9

    # EnergyWindows
    cc = sim.add_actor("DigitizerEnergyWindowsActor", "EnergyWindows")
    cc.mother = hc.mother
    cc.input_digi_collection = "Efficiency"
    cc.channels = [
        {"name": "scatter", "min": 114 * keV, "max": 126 * keV},
        {"name": "peak140", "min": 126 * keV, "max": 154.55 * keV},
        {
            "name": "spectrum",
            "min": 0 * keV,
            "max": 5000 * keV,
        },  # should be strictly equal to 'Singles'
    ]
    cc.output = hc.output

    # Modify the priority to have the correct actors order
    hc.priority = 90 # DigitizerHitsCollectionActor 
    sc.priority = 91 # DigitizerAdderActor
    bc.priority = 92 # DigitizerBlurringActor
    ea.priority = 93 # DigitizerEfficiencyActor
    cc.priority = 94 # DigitizerEnergyWindowsActor

    # ----------------------------------------
    # Run the simulation
    # ----------------------------------------
    sim.run_timing_intervals = [[0, g_runtime * sec]]

    # Start a new process to not interfere with other simulations that might have run previously.
    output = sim.run(start_new_process=True)

"""
 Main function to trigger all simulations
"""
if __name__ == "__main__":

    # Settings:
    g_activity = 60e6 #[Bq]
    g_runtime = 100 #[s]
    g_channels = 512 
    g_num_of_threads = 8
    g_visu = False
    g_distance = 20 #[cm] Distance between detector and source/phantom end
    g_collimator_type = 'lehr'
    g_isotope = "Tc99m"
    g_FWHM = 9.5 * 140.511e3/100 #[eV]
    g_plot = False

    create_simulation(g_isotope, g_collimator_type, g_activity, g_runtime, g_distance, g_channels, g_num_of_threads, g_visu, g_FWHM, g_plot)

System configuration:

I'm using a virtual server running Ubuntu 22.04.3 LTS, 2 physical processors, 16 cores, 16 threads, Intel(R) Xeon(R) Gold 6148 CPU @ 2.40 GHz, 32 GB RAM.

GATE version:

Gate version 10.0b6

See opengate_info 

    

        Attribute
        Value
    

    

        Python version
        3.11.6 (main, Oct 23 2023, 22:48:54) [GCC 11.4.0]
    

    

        Platform
        linux
    

    

        Site package
        /home/<username>/opengate_env/lib/python3.11/site-packages
    

    

        Geant4 version
        geant4-11-01-patch-01 [MT]
    

    

        Geant4 MT
        True
    

    

        Geant4 GDML
        True
    

    

        Geant4 date
        (10-February-2023)
    

    

        Geant4 data
        /home/<username>/opengate_env/lib/python3.11/site-packages/opengate_core/geant4_data
    

    

        ITK version
        5.2.1
    

    

        GATE version
        10.0b6
    

    

        GATE folder
        /home/<username>/opengate_env/lib/python3.11/site-packages/opengate
    

    

        GATE date
        2023-11-03T14:05:19 (pypi)
nkrah commented 8 months ago

Thanks for posting. And sorry for not getting back before. Can you please wait for a few days until we look at this? We are currently preparing a new beta version (beta7). Maybe that will even solve your issues.

By the way: You might consider using the current master version from the github repo via the developer install: https://opengate-python.readthedocs.io/en/latest/developer_guide.html#installation-for-developers In that way, you always have the most up-to-date version with recent fixes. The beta tags sometimes significantly lag behind the current master.

nkrah commented 8 months ago

There is now a brand new beta67 version. Just upgrade via pip and you should be good to go. Can you check if the error persists? Thanks

nkrah commented 8 months ago

From the error log you shared that G4 has actually trouble tracking the particles, potentially due to a geometry overlap, yet it states it is only a warning. Does this G4 warning coincide with cases in which your simulation crash? Or does the simulation also crash without this G4 warning?

nkrah commented 8 months ago

By the way (and you probably know this): The message "Exception: The queue is empty. The spawned process probably died." is specific to a simulation run with sim.run(start_new_process=True). If you set the keyword argument to False, you will not get this overlayed exception and only the underlying one.

franky180 commented 8 months ago

The developer install seems to be related with too many potential pitfalls I may encounter during the process of compillation. Thus, I would prefer to stick with the beta version. As suggested by @tbaudier in https://github.com/OpenGATE/opengate/issues/245#issuecomment-1748960421_ I download the latest wheels and installed them. This way, with version 10.0b6 I tought to have been up to date. I actually installed version 10.0b7 but encounter a new error message with libGL when trying to use the visualization option. The initial error still persists with this version too.

See full error message for libGL 
`libGL error: MESA-LOADER: failed to open swrast: /usr/lib/dri/swrast_dri.so: cannot open shared object file: No such file or directory (search paths /usr/lib/x86_64-linux-gnu/dri:\$${ORIGIN}/dri:/usr/lib/dri, suffix _dri)
libGL error: failed to load driver: swrast
2024-01-12 16:56:33.283 (   2.191s) [    7FA64EB33800]vtkXOpenGLRenderWindow.:651    ERR| vtkXOpenGLRenderWindow (0x5e987f0): Cannot create GLX context.  Aborting.
ERROR:root:Cannot create GLX context.  Aborting.
Aborted (core dumped)`

Concerning the G4 overlap warning, if I run the same code with less activity I don't encounter the overlap warning. I'm not able to make sense out of it.

I actually did not know that, thanks for the info with start_new_process. I only use this option to be able to run multiple simulations over night in a loop.

franky180 commented 8 months ago

Is there any way to use the default opengate/contrib/GateMaterials.db in my simulation? I wrote a custom function add_default_material_database(sim) in opengate/contrib/__init__.py in my local installation of opengate. Is there a better way to use the default GateMaterials?

Function add_default_material_database(sim) 
import pathlib
def add_default_material_database(sim):

    # Assuming your base directory is where your script or application is located
    base_directory = pathlib.Path(__file__).resolve().parent

    # Define the path to the "opengate.contrib" folder
    opengate_contrib_folder = base_directory / "opengate.contrib"

    sim.volume_manager.add_material_database(pathlib.Path(__file__).resolve().parent / "GateMaterials.db")
nkrah commented 8 months ago

Is there any way to use the default opengate/contrib/GateMaterials.db in my simulation? I wrote a custom function add_default_material_database(sim) in opengate/contrib/__init__.py in my local installation of opengate. Is there a better way to use the default GateMaterials?

Function add_default_material_database(sim)

I moved this question to a separate issue as it seems an unrelated question.

jizhang02 commented 7 months ago

Hello,

I met a similar error when I ran the alpha particle simulation. I have updated to the latest Gate version beta07. Does anyone have a solution to it? Thanks a lot.

WARNING: group: unknown groupid 1000001
terminate called after throwing an instance of 'std::out_of_range'
  what():  basic_string::substr: __pos (which is 1) > this->size() (which is 0)
/var/spool/slurm/slurmd/job124242/slurm_script: line 17: 359382 Aborted                 (core dumped) singularity run /home2/jzhang/image_gate10beta07.sif bash -c "source $HOME/set_environment.sh && python $PWD/mc_phantom_dose.py /home2/jzhang/cluster_mc/$SLURM_JOB_ID"

Here is my code:

import opengate as gate
import itk
from scipy.spatial.transform import Rotation
import numpy as np
import sys

def mc_dose_sim(path, ct_name, dose_name, source_name, total_act):

    '''
    General
    '''
    sim = gate.Simulation()
    sim.g4_verbose = True
    sim.g4_verbose_level = 0
    sim.g4_verbose_level_tracking = 1
    sim.visu = False
    sim.number_of_threads = 32
    sim.random_engine = "MersenneTwister"
    sim.random_seed = 123456
    sim.volume_manager.add_material_database(path + 'GateMaterials.db')
    m = gate.g4_units.m
    mm = gate.g4_units.mm
    um = gate.g4_units.um
    keV = gate.g4_units.keV
    MeV = gate.g4_units.MeV
    sec = gate.g4_units.second
    Bq = gate.g4_units.Bq
    kBq = 1000 * Bq
    MBq = Bq * 1e6

    '''
    Geometry
    '''
    # world
    world = sim.world
    world.size = [2 * m, 2 * m, 2 * m]
    world.material = "G4_AIR"
    # CT image
    ct = sim.add_volume('Image', 'ct')
    ct.mother = 'world'
    ct.image = path + ct_name
    ct_info = gate.image.read_image_info(ct.image)
    ct.voxel_materials = [
        [0,   4,  "G4_AIR"],
        [5,   6,  "G4_EYE_LENS_ICRP"],
        [7,   7,  "G4_SKIN_ICRP"],
        [8,  43,  "G4_BRAIN_ICRP"],
        [44, 51,  "G4_BLOOD_ICRP"],
        [52, 54,  "G4_LUNG_ICRP"],
        [55, 61,  "G4_MUSCLE_SKELETAL_ICRP"],
        [62, 63,  "G4_B-100_BONE"],
        [64, 66,  "G4_BONE_COMPACT_ICRU"],
        [67, 111, "G4_TISSUE_SOFT_ICRP"]
    ]

    '''
    Sources
    '''
    # Activity image
    source = sim.add_source("VoxelsSource", "activity")
    source.mother = ct.name
    source.image = path + source_name
    source.particle = "alpha"
    source.half_life = 857088 * sec # ac225
    source.energy.type = "mono"
    source.energy.mono = 0
    source.activity = total_act * MBq / sim.number_of_threads
    source.direction.type = "iso"
    source.position.type = "sphere"
    source.position.radius = 1 * mm
    source.position.translation = gate.image.get_translation_between_images_center(ct.image, source.image)

    '''
    Physics
    '''
    sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3"
    sim.enable_decay(True)
    sim.physics_manager.apply_cuts = True
    sim.physics_manager.set_production_cut("world", "all", 100 * um) # all means: proton, electron, positron, gamma

    '''
    Actors
    '''
    # add DoseActor
    dose = sim.add_actor('DoseActor', 'dose')
    dose.output = dose_name # dose, uncertainty, energy map
    dose.mother = 'ct'
    dose.size = ct_info.size
    dose.spacing = ct_info.spacing
    dose.img_coord_system = True 
    dose.uncertainty = True 
    dose.gray = True

    # add stat actor
    stats = sim.add_actor('SimulationStatisticsActor', 'stats')
    stats.mother = 'world'
    stats.track_types_flag = True

    '''
    Simulation
    '''
    sim.run_timing_intervals = [[0, 10 * 60 * sec]]
    sim.start()
    stats = sim.output.get_actor("stats")
    print(stats)
    print("-" * 80)

file_id     = str(1)
data_path   = '/home2/jzhang/python_code/opengate/data/'
ct_name     = 'phantom/ct/' + file_id + '.mhd'
source_name = 'phantom/act_ac/' + file_id + '.mhd'
dose_name   = str(sys.argv[1]) + '/' +file_id + '_ac.mhd'
total_act   = 133191332.4 * 1e-9 # scaled total activity MBq
mc_dose_sim(data_path, ct_name, dose_name, source_name, total_act)
nkrah commented 7 months ago

Thanks for posting. Can you share the terminal output generated by Gate/Geant4? I guess you are running the code on a cluster, but maybe you have access to the output? Otherwise, can you run and reproduce the error on a local machine where you have access to the terminal output?

jizhang02 commented 7 months ago

Thanks for posting. Can you share the terminal output generated by Gate/Geant4? I guess you are running the code on a cluster, but maybe you have access to the output? Otherwise, can you run and reproduce the error on a local machine where you have access to the terminal output? Hi, I ran on local machine and set the thread number to 1, the error disppears. Maybe the thread number can not be set too large.

Thanks a lot!

nkrah commented 7 months ago

The developer install seems to be related with too many potential pitfalls I may encounter during the process of compillation. Thus, I would prefer to stick with the beta version. As suggested by @tbaudier in #245 (comment)_ I download the latest wheels and installed them. This way, with version 10.0b6 I tought to have been up to date. I actually installed version 10.0b7 but encounter a new error message with libGL when trying to use the visualization option. The initial error still persists with this version too.

See full error message for libGL Concerning the G4 overlap warning, if I run the same code with less activity I don't encounter the overlap warning. I'm not able to make sense out of it.

I actually did not know that, thanks for the info with start_new_process. I only use this option to be able to run multiple simulations over night in a loop.

@franky180 : I was thinking: since you are having issues with the install, should we schedule a video meeting and look at a shared screen together to understand what's going on? Might be more efficient...
Drop us an email, if you want: nils.krah@creatis.insa-lyon.fr