From 4608fd05891bc9e02ce8aa1f09285e91475d9bf3 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Thu, 19 Dec 2024 11:01:00 +0100 Subject: [PATCH 01/20] first model --- garf_models/digitizers.py | 98 ++++++++++ garf_models/garf_models_helpers.py | 169 ++++++++++++++++++ garf_models/generate_garf_training_dataset.py | 66 +++++++ garf_models/readme.md | 42 +++++ garf_models/test01_helpers.py | 95 ++++++++++ .../test01_main1_intevo_source_point_ref.py | 73 ++++++++ .../test01_main2_intevo_source_point_arf.py | 68 +++++++ ...1_main3_inteov_source_point_free_flight.py | 73 ++++++++ garf_models/test01_main4_compare.py | 49 +++++ 9 files changed, 733 insertions(+) create mode 100755 garf_models/digitizers.py create mode 100755 garf_models/garf_models_helpers.py create mode 100755 garf_models/generate_garf_training_dataset.py create mode 100644 garf_models/readme.md create mode 100755 garf_models/test01_helpers.py create mode 100755 garf_models/test01_main1_intevo_source_point_ref.py create mode 100755 garf_models/test01_main2_intevo_source_point_arf.py create mode 100755 garf_models/test01_main3_inteov_source_point_free_flight.py create mode 100755 garf_models/test01_main4_compare.py diff --git a/garf_models/digitizers.py b/garf_models/digitizers.py new file mode 100755 index 0000000..a9e2077 --- /dev/null +++ b/garf_models/digitizers.py @@ -0,0 +1,98 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate as gate +from scipy.spatial.transform import Rotation + + +def add_intevo_digitizer_lu177_v3(sim, crystal_name, name, spectrum_channel=True): + + # hits + hits = sim.add_actor("DigitizerHitsCollectionActor", f"hits_{name}") + hits.attached_to = crystal_name + hits.output_filename = "" # No output + hits.attributes = [ + "PostPosition", + "TotalEnergyDeposit", + "PreStepUniqueVolumeID", + "PostStepUniqueVolumeID", + "GlobalTime", + ] + + # singles + singles = sim.add_actor("DigitizerAdderActor", f"singles_{name}") + singles.attached_to = crystal_name + singles.input_digi_collection = hits.name + # sc.policy = "EnergyWeightedCentroidPosition" + singles.policy = "EnergyWinnerPosition" + singles.output_filename = "" # No output + singles.group_volume = None + + # efficiency actor + eff = sim.add_actor("DigitizerEfficiencyActor", f"singles_{name}_eff") + eff.attached_to = crystal_name + eff.input_digi_collection = singles.name + eff.efficiency = 0.86481 # FIXME probably wrong, to evaluate + eff.efficiency = 1.0 + eff.output_filename = "" # No output + + # energy blur + keV = gate.g4_units.keV + MeV = gate.g4_units.MeV + ene_blur = sim.add_actor("DigitizerBlurringActor", f"singles_{name}_eblur") + ene_blur.output_filename = "" + ene_blur.attached_to = crystal_name + ene_blur.input_digi_collection = eff.name + ene_blur.blur_attribute = "TotalEnergyDeposit" + ene_blur.blur_method = "Linear" + ene_blur.blur_resolution = 0.13 + ene_blur.blur_reference_value = 80 * keV + ene_blur.blur_slope = -0.09 * 1 / MeV + + # spatial blurring + mm = gate.g4_units.mm + spatial_blur = sim.add_actor( + "DigitizerSpatialBlurringActor", f"singles_{name}_sblur" + ) + spatial_blur.output_filename = "" + spatial_blur.attached_to = crystal_name + spatial_blur.input_digi_collection = ene_blur.name + spatial_blur.blur_attribute = "PostPosition" + spatial_blur.blur_fwhm = 3.9 * mm + spatial_blur.keep_in_solid_limits = True + + # energy windows + singles_ene_windows = sim.add_actor( + "DigitizerEnergyWindowsActor", f"{name}_energy_window" + ) + channels = [ + {"name": f"spectrum", "min": 3 * keV, "max": 515 * keV}, + {"name": f"scatter1_{name}", "min": 96 * keV, "max": 104 * keV}, + {"name": f"peak113_{name}", "min": 104.52 * keV, "max": 121.48 * keV}, + {"name": f"scatter2_{name}", "min": 122.48 * keV, "max": 133.12 * keV}, + {"name": f"scatter3_{name}", "min": 176.46 * keV, "max": 191.36 * keV}, + {"name": f"peak208_{name}", "min": 192.4 * keV, "max": 223.6 * keV}, + {"name": f"scatter4_{name}", "min": 224.64 * keV, "max": 243.3 * keV}, + ] + if not spectrum_channel: + channels.pop(0) + singles_ene_windows.attached_to = crystal_name + singles_ene_windows.input_digi_collection = spatial_blur.name + singles_ene_windows.channels = channels + singles_ene_windows.output_filename = "" # No output + + # projection + proj = sim.add_actor("DigitizerProjectionActor", f"{name}_projection") + proj.attached_to = crystal_name + proj.input_digi_collections = [x["name"] for x in channels] + proj.spacing = [4.7951998710632 * mm / 2, 4.7951998710632 * mm / 2] + proj.size = [128 * 2, 128 * 2] + proj.output_filename = "projection.mhd" + proj.origin_as_image_center = True + + # plane orientation + proj.detector_orientation_matrix = Rotation.from_euler( + "yx", (90, 90), degrees=True + ).as_matrix() + + return proj diff --git a/garf_models/garf_models_helpers.py b/garf_models/garf_models_helpers.py new file mode 100755 index 0000000..d07bda6 --- /dev/null +++ b/garf_models/garf_models_helpers.py @@ -0,0 +1,169 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate.contrib.spect.ge_discovery_nm670 as nm670 +import opengate.contrib.spect.siemens_intevo as intevo +from opengate.sources.base import get_rad_gamma_spectrum +from opengate.exception import fatal +from digitizers import * + +all_spects = ("intevo", "nm670") + +all_collimators = { + "intevo": {"tc99m": "lehr", "lu177": "melp", "in111": "melp", "i131": "he"}, + "nm670": {"tc99m": "lehr", "lu177": "megp", "in111": "megp", "i131": "hegp"}, +} + +all_digitizers = { + "intevo": { + "v1": { + "tc99m": intevo.add_digitizer_tc99m, + "lu177": intevo.add_digitizer_lu177, + }, + "v2": {"tc99m": intevo.add_digitizer_tc99m_v2}, + "v3": {"lu177": add_intevo_digitizer_lu177_v3}, + }, + "nm670": { + "v1": { + "tc99m": nm670.add_digitizer_tc99m, + "lu177": nm670.add_digitizer_lu177, + }, + "v2": {"tc99m": nm670.add_digitizer_tc99m_v2}, + }, +} + + +def init_default_garf_simulation(sim, visu, threads): + # main options + sim.visu = visu + sim.visu_type = "qt" + sim.number_of_threads = threads + sim.progress_bar = True + sim.store_json_archive = True + sim.store_input_files = False + sim.json_archive_filename = "simu.json" + + # units + mm = gate.g4_units.mm + m = gate.g4_units.m + + # world + world = sim.world + world.size = [2 * m, 2 * m, 2 * m] + world.material = "G4_AIR" + + # physics + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option4" + sim.physics_manager.global_production_cuts.all = 1 * mm + + # add stat actor + stats = sim.add_actor("SimulationStatisticsActor", "stats") + stats.output_filename = "stats.txt" + return stats + + +def get_collimator_from_rad(spect, rad): + spect = spect.lower() + rad = rad.lower() + if spect not in all_collimators: + fatal(f'Unknown spect system "{spect}", known are: {all_collimators.keys()} ') + if rad not in all_collimators[spect]: + fatal( + f'Unknown radionuclide "{rad}", known are: {all_collimators[spect].keys()} ' + ) + return all_collimators[spect][rad] + + +def add_spect_imaging_device(sim, spect, rad, crystal_size): + if spect not in all_spects: + fatal(f'Unknown spect system "{spect}", known are: {all_spects} ') + colli_type = get_collimator_from_rad(spect, rad) + mm = gate.g4_units.mm + head, colli, crystal = None, None, None + if spect == "intevo": + head, colli, crystal = intevo.add_spect_head( + sim, "spect", collimator_type=colli_type, debug=sim.visu == True + ) + intevo.set_head_orientation(head, colli_type, radius=0 * mm) + if spect == "nm670": + head, colli, crystal = nm670.add_spect_head( + sim, + "spect", + collimator_type=colli_type, + debug=sim.visu == True, + crystal_size=crystal_size, + ) + + print("Head position", head.translation) + return head, colli, crystal, colli_type + + +def add_digitizer(sim, spect, rad, digitizer, crystal): + if spect not in all_digitizers: + fatal( + f'Unknown digitizer spect system "{spect}", known are: {all_digitizers.keys()} ' + ) + digitizers = all_digitizers[spect] + if digitizer not in digitizers: + fatal(f'Unknown digitizer "{digitizer}", known are: {digitizers.keys()} ') + digitizers = digitizers[digitizer] + if rad not in digitizers: + fatal( + f'Unknown digitizer "{digitizer}" for radionuclide "{rad}", known are: {digitizers.keys()} ' + ) + + # create the digitizer + f = digitizers[rad] + + # digitizer + f(sim, crystal, f"digitizer") + + # special case for spectrum channel (to remove) + proj = sim.get_actor(f"digitizer_projection") + proj.write_to_disk = False + ew = sim.get_actor(f"digitizer_energy_window") + channels = ew.channels + c = [] + for channel in channels: + if channel.name != "spectrum": + c.append(channel) + ew.channels = c + proj.input_digi_collections = [c["name"] for c in ew.channels] + + return ew + + +def add_arf(sim, spect, head, colli_type, ew, rr): + detector_plane = None + arf = None + if spect not in all_spects: + fatal(f'Unknown spect system "{spect}", known are: {all_spects} ') + if spect == "nm670": + detector_plane, arf = nm670.add_actor_for_arf_training_dataset( + sim, head, colli_type, ew, rr=rr + ) + if spect == "intevo": + detector_plane, arf = intevo.add_actor_for_arf_training_dataset( + sim, colli_type, ew, rr=rr + ) + print("Plane position", detector_plane.translation) + return detector_plane, arf + + +def add_source(sim, spect, rad, activity, detector_plane): + if spect not in all_spects: + fatal(f'Unknown spect system "{spect}", known are: {all_spects} ') + MeV = gate.g4_units.MeV + max_energy = get_rad_gamma_spectrum(rad).energies[-1] * 1.05 + print(f"Max energy for {rad} is {max_energy} MeV") + source = None + if spect == "nm670": + source = nm670.add_source_for_arf_training_dataset( + sim, "src", activity, detector_plane, 0.01 * MeV, max_energy + ) + if spect == "intevo": + source = intevo.add_source_for_arf_training_dataset( + sim, "src", activity, detector_plane, 0.01 * MeV, max_energy + ) + print("Source position", source.position.translation) + return source diff --git a/garf_models/generate_garf_training_dataset.py b/garf_models/generate_garf_training_dataset.py new file mode 100755 index 0000000..f2f79b6 --- /dev/null +++ b/garf_models/generate_garf_training_dataset.py @@ -0,0 +1,66 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +from pathlib import Path +import click +from garf_models_helpers import * + +CONTEXT_SETTINGS = dict(help_option_names=["-h", "--help"]) + + +@click.command(context_settings=CONTEXT_SETTINGS) +@click.option("--spect", "-s", default="intevo", help="SPECT system: intevo, nm670") +@click.option("--rad", "-r", required=True, help="Radionuclide 177lu, 99tcm etc") +@click.option("--digitizer", "-d", required=True, help="Digitizer version : v1 v2 etc") +@click.option("--n", "-n", default=5e8, help="Number of particles") +@click.option("--rr", default=50, help="Russian Roulette for GARF") +@click.option("--visu", is_flag=True, default=False, help="visu only") +@click.option("--threads", "-t", default=4, help="number of threads") +@click.option( + "--crystal_size", "-c", default="5/8", help="For nm670 crystal size 3/8 or 5/8" +) +def go(spect, rad, digitizer, n, rr, visu, threads, crystal_size): + + # set the default simulation param + spect_type = spect + sim = gate.Simulation() + arf_basename = f"{spect_type}_{rad}_{digitizer}" + stats = init_default_garf_simulation(sim, visu, threads) + sim.output_dir = Path("output") / arf_basename + + # set the SPECT system + print(f"SPECT type is : {spect_type}") + head, colli, crystal, colli_type = add_spect_imaging_device( + sim, spect_type, rad, crystal_size + ) + print(f"collimator type is : {colli_type}") + + # set the digitizer and the arf + ew = add_digitizer(sim, spect_type, rad, digitizer, crystal) + detector_plane, arf = add_arf(sim, spect_type, head, colli_type, ew, rr) + + # visu option + if visu: + n = 1e2 + sim.number_of_threads = 1 + sim.output_dir = Path("output") / f"visu" + + # set the source + Bq = gate.g4_units.Bq + activity = n * Bq / sim.number_of_threads + add_source(sim, spect_type, rad, activity, detector_plane) + + # go + sim.run() + + # print results at the end + print(stats) + + # print cmd line to train GARF + pth_json = Path("pth") / "train_arf_v034.json" + pth = Path("pth") / f"{arf_basename}.pth" + print(f"garf_train {pth_json} {arf.get_output_path()} {pth}") + + +# -------------------------------------------------------------------------- +if __name__ == "__main__": + go() diff --git a/garf_models/readme.md b/garf_models/readme.md new file mode 100644 index 0000000..1ff3a11 --- /dev/null +++ b/garf_models/readme.md @@ -0,0 +1,42 @@ + + + +## test01 + +- Intevo +- Lu177, melp +- 2 point-sources + AA + + +## garf generation + + ./generate_garf_training_dataset.py -s intevo -r tc99m -d v1 + ./generate_garf_training_dataset.py -s intevo -r lu177 -d v1 + ./generate_garf_training_dataset.py -s intevo -r tc99m -d v2 + + # WIP expected 12 min 200 MB + ./generate_garf_training_dataset.py -s intevo -r lu177 -d v3 -n 2e9 + + ./generate_garf_training_dataset.py -s nm670 -r tc99m -d v1 + ./generate_garf_training_dataset.py -s nm670 -r tc99m -d v2 + ./generate_garf_training_dataset.py -s nm670 -r lu177 -d v1 + + +## Ideal UI ? + +options +- device: intevo, nm670 +- collimator: lehr, etc +- rad: 177lu, Tc99m etc +- digit + windowing: ?????? +- physics, cuts +- rr=50 ? +- nb particle, root size +- Output: stat + root + simu.json + pth + +Naming: +- intevo_melp_177lu_digit_2_v034.pth +- nm670_3p8_lehr_99tcm_digit_33_v034.pth +- nm670_5p8_lehr_99tcm_digit_33_v034.pth + + generate_training_dataset -s intevo -r 177lu --digit x diff --git a/garf_models/test01_helpers.py b/garf_models/test01_helpers.py new file mode 100755 index 0000000..6de0b55 --- /dev/null +++ b/garf_models/test01_helpers.py @@ -0,0 +1,95 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate.contrib.phantoms.nemaiec as nemaiec +from opengate.sources.base import set_source_rad_energy_spectrum +from digitizers import * +import numpy as np + +def init_sim(sim): + m = gate.g4_units.m + mm = gate.g4_units.mm + # world + world = sim.world + world.size = [2 * m, 2 * m, 2 * m] + world.material = "G4_AIR" + # physics + sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3" + sim.physics_manager.global_production_cuts.all = 1 * mm + # add stat actor + stats = sim.add_actor("SimulationStatisticsActor", "stats") + stats.output_filename = f"stats.txt" + return stats + + + +def add_vox_iec(sim, spacing, data_path): + mm = gate.g4_units.mm + # voxelize_iec_phantom -o data/iec_5mm.mha --spacing 5 + print(f"Phantom: IEC voxelized {spacing}mm") + mhd_filename = data_path / f"iec_{spacing}mm.mha" + labels_filename = data_path / f"iec_{spacing}mm.json" + iec = sim.add_volume("Image", "iec") + iec.image = mhd_filename + nemaiec.create_material(sim) + iec.set_materials_from_voxelisation(labels_filename) + sim.physics_manager.set_production_cut(iec.name, "all", 2 * mm) + return iec + + +def add_vox_source(sim, rad, activity, data_path): + mm = gate.g4_units.mm + Bq = gate.g4_units.Bq + # voxelize_iec_phantom -o data/iec_1mm.mhd --spacing 1 --output_source data/iec_1mm_activity.mhd -a 1 2 3 4 5 6 + iec_source_filename = data_path / "iec_1mm_activity.mhd" + source = sim.add_source("VoxelSource", "src") + source.image = iec_source_filename + source.position.translation = [0, 35 * mm, 0] + set_source_rad_energy_spectrum(source, rad) + source.particle = "gamma" + _, volumes = nemaiec.get_default_sphere_centers_and_volumes() + print(f"Volumes are {volumes}") + source.activity = activity * np.array(volumes).sum() + print(f"Total activity is {source.activity / Bq}") + return source + +def add_source_point_test(sim, rad, planes, n, angle_tolerance, head_type='arf'): + keV = gate.g4_units.keV + source = sim.add_source("GenericSource", "source_point") + set_source_rad_energy_spectrum(source, rad) + source.particle = "gamma" + source.direction.type = "iso" + source.direction.acceptance_angle.volumes = [p.name for p in planes] + source.direction.acceptance_angle.skip_policy = "SkipEvents" + source.direction.acceptance_angle.intersection_flag = True + source.direction.acceptance_angle.normal_flag = True + source.direction.acceptance_angle.normal_vector = [0, 0, -1] + if head_type == 'spect': + source.direction.acceptance_angle.normal_vector = [1, 0, 0] + source.direction.acceptance_angle.normal_tolerance = angle_tolerance + source.n = n + source.attached_to = 'b' + + cm = gate.g4_units.cm + + b = sim.add_volume('Box', 'b') + b.translation = [3*cm, 10*cm, 2*cm] + b.translation = [3*cm, 0*cm, 2*cm] + b.material = 'G4_AIR' + + source = sim.add_source("GenericSource", "source_point2") + set_source_rad_energy_spectrum(source, rad) + source.particle = "gamma" + source.direction.type = "iso" + source.direction.acceptance_angle.volumes = [p.name for p in planes] + source.direction.acceptance_angle.skip_policy = "SkipEvents" + source.direction.acceptance_angle.intersection_flag = True + source.direction.acceptance_angle.normal_flag = True + source.direction.acceptance_angle.normal_vector = [0, 0, -1] + if head_type == 'spect': + source.direction.acceptance_angle.normal_vector = [1, 0, 0] + source.direction.acceptance_angle.normal_tolerance = angle_tolerance + source.n = n/2 + + + return source \ No newline at end of file diff --git a/garf_models/test01_main1_intevo_source_point_ref.py b/garf_models/test01_main1_intevo_source_point_ref.py new file mode 100755 index 0000000..5ae1e39 --- /dev/null +++ b/garf_models/test01_main1_intevo_source_point_ref.py @@ -0,0 +1,73 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate.contrib.spect.siemens_intevo as intevo +from pathlib import Path +from digitizers import * +from test01_helpers import init_sim, add_source_point_test + + +def go(): + + # create the simulation + sim = gate.Simulation() + + # main options + #sim.visu = True + sim.visu_type = "qt" + sim.random_seed = "auto" + sim.number_of_threads = 4 + sim.progress_bar = True + sim.output_dir = Path("test01") / f"reference" + + # units + cm = gate.g4_units.cm + deg = gate.g4_units.deg + + # options + radius = 28 * cm + rad = "lu177" + colli_type = "melp" + activity = 4e5 + angle_tolerance = 10 * deg + + # visu + if sim.visu: + sim.number_of_threads = 1 + activity = 100 + + # world etc + stats = init_sim(sim) + + # set the spect head + head1, colli1, crystal1 = intevo.add_spect_head( + sim, "spect1", collimator_type=colli_type, debug=sim.visu == True + ) + proj1 = add_intevo_digitizer_lu177_v3( + sim, crystal1, f"digitizer1", spectrum_channel=False + ) + proj1.output_filename = "projection_1.mhd" + + head2, colli2, crystal2 = intevo.add_spect_head( + sim, "spect2", collimator_type=colli_type, debug=sim.visu == True + ) + proj2 = add_intevo_digitizer_lu177_v3( + sim, crystal2, f"digitizer2", spectrum_channel=False + ) + proj2.output_filename = "projection_2.mhd" + + # rotate + intevo.rotate_gantry(head1, radius, 0, initial_rotation='spect') + intevo.rotate_gantry(head2, radius, 180, initial_rotation='spect') + + # source for test + add_source_point_test(sim, rad, (head1, head2), activity, angle_tolerance, head_type='spect') + + # go + sim.run() + print(stats) + + +# -------------------------------------------------------------------------- +if __name__ == "__main__": + go() diff --git a/garf_models/test01_main2_intevo_source_point_arf.py b/garf_models/test01_main2_intevo_source_point_arf.py new file mode 100755 index 0000000..19ec0cd --- /dev/null +++ b/garf_models/test01_main2_intevo_source_point_arf.py @@ -0,0 +1,68 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate.contrib.spect.siemens_intevo as intevo +from pathlib import Path +from digitizers import * +from test01_helpers import init_sim, add_source_point_test + + +def go(): + + # create the simulation + sim = gate.Simulation() + + # main options + #sim.visu = True + sim.visu_type = "qt" + sim.random_seed = "auto" + sim.number_of_threads = 4 + sim.progress_bar = True + sim.output_dir = Path("test01") / f"arf" + + # units + mm = gate.g4_units.mm + cm = gate.g4_units.cm + deg = gate.g4_units.deg + + # options + radius = 28 * cm + rad = "lu177" + colli_type = "melp" + activity = 4e5 + angle_tolerance = 10 * deg + + # visu + if sim.visu: + sim.number_of_threads = 1 + activity = 100 + + # world etc + stats = init_sim(sim) + + # set the two spect heads + spacing = [4.7951998710632 * mm / 2, 4.7951998710632 * mm / 2] + size = [128 * 2, 128 * 2] + pth = Path("pth") / "intevo_lu177_v3.pth" + det_plane1, arf1 = intevo.add_arf_detector( + sim, radius, 0, size, spacing, colli_type, "detector", 1, pth + ) + det_plane2, arf2 = intevo.add_arf_detector( + sim, radius, 180, size, spacing, colli_type, "detector", 2, pth + ) + # compute the gantry rotations + intevo.rotate_gantry(det_plane1, radius, 0, initial_rotation="arf") + intevo.rotate_gantry(det_plane2, radius, 180, initial_rotation="arf") + det_planes = [det_plane1, det_plane2] + + # source for test + add_source_point_test(sim, rad, det_planes, activity, angle_tolerance) + + # go + sim.run() + print(stats) + + +# -------------------------------------------------------------------------- +if __name__ == "__main__": + go() diff --git a/garf_models/test01_main3_inteov_source_point_free_flight.py b/garf_models/test01_main3_inteov_source_point_free_flight.py new file mode 100755 index 0000000..4d4cf35 --- /dev/null +++ b/garf_models/test01_main3_inteov_source_point_free_flight.py @@ -0,0 +1,73 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate.contrib.spect.siemens_intevo as intevo +from pathlib import Path +from digitizers import * +from test01_helpers import init_sim, add_source_point_test + + +def go(): + + # create the simulation + sim = gate.Simulation() + + # main options + #sim.visu = True + sim.visu_type = "qt" + sim.random_seed = "auto" + sim.number_of_threads = 4 + sim.progress_bar = True + sim.output_dir = Path("test01") / f"free_flight" + + # units + mm = gate.g4_units.mm + cm = gate.g4_units.cm + deg = gate.g4_units.deg + + # options + radius = 28 * cm + rad = "lu177" + colli_type = "melp" + activity = 4e5 + angle_tolerance = 10 * deg + + # visu + if sim.visu: + sim.number_of_threads = 1 + activity = 100 + + # world etc + stats = init_sim(sim) + + # set the two spect heads + spacing = [4.7951998710632 * mm / 2, 4.7951998710632 * mm / 2] + size = [128 * 2, 128 * 2] + pth = Path("pth") / "intevo_lu177_v3.pth" + det_plane1, arf1 = intevo.add_arf_detector( + sim, radius, 0, size, spacing, colli_type, "detector", 1, pth + ) + det_plane2, arf2 = intevo.add_arf_detector( + sim, radius, 180, size, spacing, colli_type, "detector", 2, pth + ) + # compute the gantry rotations + intevo.rotate_gantry(det_plane1, radius, 0, initial_rotation="arf") + intevo.rotate_gantry(det_plane2, radius, 180, initial_rotation="arf") + det_planes = [det_plane1, det_plane2] + + # source for test + add_source_point_test(sim, rad, det_planes, activity, angle_tolerance) + + # FF + ff = sim.add_actor("FreeFlightActor", "ff") + ff.attached_to = "world" + ff.particles = "gamma" + + # go + sim.run() + print(stats) + + +# -------------------------------------------------------------------------- +if __name__ == "__main__": + go() diff --git a/garf_models/test01_main4_compare.py b/garf_models/test01_main4_compare.py new file mode 100755 index 0000000..2d2a5a8 --- /dev/null +++ b/garf_models/test01_main4_compare.py @@ -0,0 +1,49 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate.contrib.spect.siemens_intevo as intevo +from pathlib import Path +from digitizers import * +from test01_helpers import init_sim, add_source_point_test +from opengate.tests import utility + + +def go(): + + folder_base = Path('test01') + folder_ref = folder_base / "reference" + folder_arf = folder_base / "arf" + folder_ff = folder_base / "free_flight" + + stats_ref = utility.read_stat_file(folder_ref / "stats.txt") + + is_ok = True + + im_name = "projection_1.mhd" + is_ok = utility.assert_images( + folder_ref / im_name, + folder_arf / im_name, + stats_ref, + tolerance=38, + ignore_value_data1=0, + sum_tolerance=8, + axis="x", + ) and is_ok + + im_name = "projection_2.mhd" + is_ok = utility.assert_images( + folder_ref / im_name, + folder_arf / im_name, + stats_ref, + tolerance=38, + ignore_value_data1=0, + sum_tolerance=8, + axis="x", + ) and is_ok + utility.test_ok(is_ok) + + + +# -------------------------------------------------------------------------- +if __name__ == "__main__": + go() From 1c2484a24454c0d12c68dd3af8c3661ed7c00b8c Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Thu, 19 Dec 2024 11:39:02 +0100 Subject: [PATCH 02/20] compare --- garf_models/test01_helpers.py | 19 +++-- .../test01_main1_intevo_source_point_ref.py | 6 +- .../test01_main2_intevo_source_point_arf.py | 6 +- ...1_main3_inteov_source_point_free_flight.py | 73 ------------------- 4 files changed, 17 insertions(+), 87 deletions(-) delete mode 100755 garf_models/test01_main3_inteov_source_point_free_flight.py diff --git a/garf_models/test01_helpers.py b/garf_models/test01_helpers.py index 6de0b55..916b3a1 100755 --- a/garf_models/test01_helpers.py +++ b/garf_models/test01_helpers.py @@ -53,12 +53,14 @@ def add_vox_source(sim, rad, activity, data_path): print(f"Total activity is {source.activity / Bq}") return source -def add_source_point_test(sim, rad, planes, n, angle_tolerance, head_type='arf'): - keV = gate.g4_units.keV +def add_source_point_test(sim, rad, planes, activity, angle_tolerance, head_type='arf'): + cm = gate.g4_units.cm source = sim.add_source("GenericSource", "source_point") set_source_rad_energy_spectrum(source, rad) source.particle = "gamma" source.direction.type = "iso" + source.position.type = "sphere" + source.position.radius = 1*cm source.direction.acceptance_angle.volumes = [p.name for p in planes] source.direction.acceptance_angle.skip_policy = "SkipEvents" source.direction.acceptance_angle.intersection_flag = True @@ -67,13 +69,11 @@ def add_source_point_test(sim, rad, planes, n, angle_tolerance, head_type='arf') if head_type == 'spect': source.direction.acceptance_angle.normal_vector = [1, 0, 0] source.direction.acceptance_angle.normal_tolerance = angle_tolerance - source.n = n + source.activity = activity source.attached_to = 'b' - cm = gate.g4_units.cm - b = sim.add_volume('Box', 'b') - b.translation = [3*cm, 10*cm, 2*cm] + #b.translation = [3*cm, 10*cm, 2*cm] b.translation = [3*cm, 0*cm, 2*cm] b.material = 'G4_AIR' @@ -81,6 +81,8 @@ def add_source_point_test(sim, rad, planes, n, angle_tolerance, head_type='arf') set_source_rad_energy_spectrum(source, rad) source.particle = "gamma" source.direction.type = "iso" + source.position.type = "sphere" + source.position.radius = 1*cm source.direction.acceptance_angle.volumes = [p.name for p in planes] source.direction.acceptance_angle.skip_policy = "SkipEvents" source.direction.acceptance_angle.intersection_flag = True @@ -89,7 +91,4 @@ def add_source_point_test(sim, rad, planes, n, angle_tolerance, head_type='arf') if head_type == 'spect': source.direction.acceptance_angle.normal_vector = [1, 0, 0] source.direction.acceptance_angle.normal_tolerance = angle_tolerance - source.n = n/2 - - - return source \ No newline at end of file + source.activity = activity/2 \ No newline at end of file diff --git a/garf_models/test01_main1_intevo_source_point_ref.py b/garf_models/test01_main1_intevo_source_point_ref.py index 5ae1e39..3112725 100755 --- a/garf_models/test01_main1_intevo_source_point_ref.py +++ b/garf_models/test01_main1_intevo_source_point_ref.py @@ -23,18 +23,19 @@ def go(): # units cm = gate.g4_units.cm deg = gate.g4_units.deg + Bq = gate.g4_units.Bq # options radius = 28 * cm rad = "lu177" colli_type = "melp" - activity = 4e5 + activity = 1e8 * Bq angle_tolerance = 10 * deg # visu if sim.visu: sim.number_of_threads = 1 - activity = 100 + activity = 1000 * Bq # world etc stats = init_sim(sim) @@ -57,6 +58,7 @@ def go(): proj2.output_filename = "projection_2.mhd" # rotate + nb_angle = 3 intevo.rotate_gantry(head1, radius, 0, initial_rotation='spect') intevo.rotate_gantry(head2, radius, 180, initial_rotation='spect') diff --git a/garf_models/test01_main2_intevo_source_point_arf.py b/garf_models/test01_main2_intevo_source_point_arf.py index 19ec0cd..bce0b4c 100755 --- a/garf_models/test01_main2_intevo_source_point_arf.py +++ b/garf_models/test01_main2_intevo_source_point_arf.py @@ -24,18 +24,19 @@ def go(): mm = gate.g4_units.mm cm = gate.g4_units.cm deg = gate.g4_units.deg + Bq = gate.g4_units.Bq # options radius = 28 * cm rad = "lu177" colli_type = "melp" - activity = 4e5 + activity = 1e8* Bq angle_tolerance = 10 * deg # visu if sim.visu: sim.number_of_threads = 1 - activity = 100 + activity = 1000* Bq # world etc stats = init_sim(sim) @@ -51,6 +52,7 @@ def go(): sim, radius, 180, size, spacing, colli_type, "detector", 2, pth ) # compute the gantry rotations + nb_angle = 3 intevo.rotate_gantry(det_plane1, radius, 0, initial_rotation="arf") intevo.rotate_gantry(det_plane2, radius, 180, initial_rotation="arf") det_planes = [det_plane1, det_plane2] diff --git a/garf_models/test01_main3_inteov_source_point_free_flight.py b/garf_models/test01_main3_inteov_source_point_free_flight.py deleted file mode 100755 index 4d4cf35..0000000 --- a/garf_models/test01_main3_inteov_source_point_free_flight.py +++ /dev/null @@ -1,73 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -import opengate.contrib.spect.siemens_intevo as intevo -from pathlib import Path -from digitizers import * -from test01_helpers import init_sim, add_source_point_test - - -def go(): - - # create the simulation - sim = gate.Simulation() - - # main options - #sim.visu = True - sim.visu_type = "qt" - sim.random_seed = "auto" - sim.number_of_threads = 4 - sim.progress_bar = True - sim.output_dir = Path("test01") / f"free_flight" - - # units - mm = gate.g4_units.mm - cm = gate.g4_units.cm - deg = gate.g4_units.deg - - # options - radius = 28 * cm - rad = "lu177" - colli_type = "melp" - activity = 4e5 - angle_tolerance = 10 * deg - - # visu - if sim.visu: - sim.number_of_threads = 1 - activity = 100 - - # world etc - stats = init_sim(sim) - - # set the two spect heads - spacing = [4.7951998710632 * mm / 2, 4.7951998710632 * mm / 2] - size = [128 * 2, 128 * 2] - pth = Path("pth") / "intevo_lu177_v3.pth" - det_plane1, arf1 = intevo.add_arf_detector( - sim, radius, 0, size, spacing, colli_type, "detector", 1, pth - ) - det_plane2, arf2 = intevo.add_arf_detector( - sim, radius, 180, size, spacing, colli_type, "detector", 2, pth - ) - # compute the gantry rotations - intevo.rotate_gantry(det_plane1, radius, 0, initial_rotation="arf") - intevo.rotate_gantry(det_plane2, radius, 180, initial_rotation="arf") - det_planes = [det_plane1, det_plane2] - - # source for test - add_source_point_test(sim, rad, det_planes, activity, angle_tolerance) - - # FF - ff = sim.add_actor("FreeFlightActor", "ff") - ff.attached_to = "world" - ff.particles = "gamma" - - # go - sim.run() - print(stats) - - -# -------------------------------------------------------------------------- -if __name__ == "__main__": - go() From 500af1334d81ea1a0b08c56b4a1f02ad73a4738c Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Fri, 20 Dec 2024 08:20:43 +0100 Subject: [PATCH 03/20] detector orientation --- garf_models/garf_models_helpers.py | 46 ++++++++++- garf_models/generate_garf_training_dataset.py | 4 +- garf_models/old/test01_main4_compare.py | 54 +++++++++++++ garf_models/pth/train_arf_v034_test.json | 33 ++++++++ garf_models/readme.md | 26 +------ garf_models/test01_helpers.py | 28 +++---- .../test01_main1_intevo_source_point_ref.py | 19 +++-- .../test01_main2_intevo_source_point_arf.py | 23 +++--- ...1_main3_intevo_source_point_free_flight.py | 78 +++++++++++++++++++ garf_models/test01_main4_compare.py | 77 ++++++++++++------ 10 files changed, 307 insertions(+), 81 deletions(-) create mode 100755 garf_models/old/test01_main4_compare.py create mode 100644 garf_models/pth/train_arf_v034_test.json create mode 100755 garf_models/test01_main3_intevo_source_point_free_flight.py diff --git a/garf_models/garf_models_helpers.py b/garf_models/garf_models_helpers.py index d07bda6..95663cf 100755 --- a/garf_models/garf_models_helpers.py +++ b/garf_models/garf_models_helpers.py @@ -74,6 +74,30 @@ def get_collimator_from_rad(spect, rad): return all_collimators[spect][rad] +def add_spect_imaging_device_OLD_TO_REMOVE(sim, spect, rad, crystal_size): + if spect not in all_spects: + fatal(f'Unknown spect system "{spect}", known are: {all_spects} ') + colli_type = get_collimator_from_rad(spect, rad) + mm = gate.g4_units.mm + head, colli, crystal = None, None, None + if spect == "intevo": + head, colli, crystal = intevo.add_spect_head( + sim, "spect", collimator_type=colli_type, debug=sim.visu == True + ) + intevo.set_head_orientation_OLD_TO_REMOVE(head, colli_type, radius=0 * mm) + if spect == "nm670": + head, colli, crystal = nm670.add_spect_head( + sim, + "spect", + collimator_type=colli_type, + debug=sim.visu == True, + crystal_size=crystal_size, + ) + + print("Head position", head.translation) + return head, colli, crystal, colli_type + + def add_spect_imaging_device(sim, spect, rad, crystal_size): if spect not in all_spects: fatal(f'Unknown spect system "{spect}", known are: {all_spects} ') @@ -84,7 +108,8 @@ def add_spect_imaging_device(sim, spect, rad, crystal_size): head, colli, crystal = intevo.add_spect_head( sim, "spect", collimator_type=colli_type, debug=sim.visu == True ) - intevo.set_head_orientation(head, colli_type, radius=0 * mm) + # intevo.set_head_orientation_OLD_TO_REMOVE(head, colli_type, radius=0 * mm) + intevo.rotate_gantry(head, radius=0, start_angle_deg=0) if spect == "nm670": head, colli, crystal = nm670.add_spect_head( sim, @@ -133,7 +158,24 @@ def add_digitizer(sim, spect, rad, digitizer, crystal): return ew -def add_arf(sim, spect, head, colli_type, ew, rr): +def add_arf_OLD_TO_REMOVE(sim, spect, head, colli_type, ew, rr): + detector_plane = None + arf = None + if spect not in all_spects: + fatal(f'Unknown spect system "{spect}", known are: {all_spects} ') + if spect == "nm670": + detector_plane, arf = nm670.add_actor_for_arf_training_dataset( + sim, head, colli_type, ew, rr=rr + ) + if spect == "intevo": + detector_plane, arf = intevo.add_actor_for_arf_training_dataset( + sim, colli_type, ew, rr=rr + ) + print("Plane position", detector_plane.translation) + return detector_plane, arf + + +def add_arf_training(sim, spect, head, colli_type, ew, rr): detector_plane = None arf = None if spect not in all_spects: diff --git a/garf_models/generate_garf_training_dataset.py b/garf_models/generate_garf_training_dataset.py index f2f79b6..fbc6d1d 100755 --- a/garf_models/generate_garf_training_dataset.py +++ b/garf_models/generate_garf_training_dataset.py @@ -36,7 +36,9 @@ def go(spect, rad, digitizer, n, rr, visu, threads, crystal_size): # set the digitizer and the arf ew = add_digitizer(sim, spect_type, rad, digitizer, crystal) - detector_plane, arf = add_arf(sim, spect_type, head, colli_type, ew, rr) + + # set the arf plane for training + detector_plane, arf = add_arf_training(sim, spect_type, head, colli_type, ew, rr) # visu option if visu: diff --git a/garf_models/old/test01_main4_compare.py b/garf_models/old/test01_main4_compare.py new file mode 100755 index 0000000..65a997d --- /dev/null +++ b/garf_models/old/test01_main4_compare.py @@ -0,0 +1,54 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate.contrib.spect.siemens_intevo as intevo +from pathlib import Path +from digitizers import * +from test01_helpers import init_sim, add_source_point_test +from opengate.tests import utility + + +def go(): + + folder_base = Path("../test01") + folder_ref = folder_base / "reference" + folder_arf = folder_base / "arf" + folder_ff = folder_base / "free_flight" + + stats_ref = utility.read_stat_file(folder_ref / "stats.txt") + + is_ok = True + + im_name = "projection_1.mhd" + is_ok = ( + utility.assert_images( + folder_ref / im_name, + folder_arf / im_name, + stats_ref, + tolerance=38, + ignore_value_data1=0, + sum_tolerance=8, + axis="x", + ) + and is_ok + ) + + im_name = "projection_2.mhd" + is_ok = ( + utility.assert_images( + folder_ref / im_name, + folder_arf / im_name, + stats_ref, + tolerance=38, + ignore_value_data1=0, + sum_tolerance=8, + axis="x", + ) + and is_ok + ) + utility.test_ok(is_ok) + + +# -------------------------------------------------------------------------- +if __name__ == "__main__": + go() diff --git a/garf_models/pth/train_arf_v034_test.json b/garf_models/pth/train_arf_v034_test.json new file mode 100644 index 0000000..0b975ef --- /dev/null +++ b/garf_models/pth/train_arf_v034_test.json @@ -0,0 +1,33 @@ +{ + "#": "Seed for random generator", + "seed": 123654, + + "#": "Russian Roulette param", + "RR": 50, + + "#": "Number of neurons by layer", + "H": 200, + + "#": "Number of layer", + "L": 3, + + "#": "Learning rate", + "learning_rate": 0.00001, + + "#": "Max nb of epoch (iteration)", + "epoch_max": 100, + "#epoch_max": 4, + + "#": "Nb of epoch without progress to stop earlier", + "early_stopping": 5, + + "#": "Nb of samples by batch", + "batch_size": 2000, + + "#": "Nb of batch per epoch NOT USED !!!!", + "batch_per_epoch": 200, + + "#": "Every epoch: store the current model (in epoch)", + "epoch_store_every": 50000 + +} diff --git a/garf_models/readme.md b/garf_models/readme.md index 1ff3a11..7ac449c 100644 --- a/garf_models/readme.md +++ b/garf_models/readme.md @@ -5,38 +5,20 @@ - Intevo - Lu177, melp +- no phantom (in Air) - 2 point-sources + AA ## garf generation + # osx = 22 min, 346 MB + ./generate_garf_training_dataset.py -s intevo -r lu177 -d v3 -n 2e9 -t 4 + ./generate_garf_training_dataset.py -s intevo -r tc99m -d v1 ./generate_garf_training_dataset.py -s intevo -r lu177 -d v1 ./generate_garf_training_dataset.py -s intevo -r tc99m -d v2 - # WIP expected 12 min 200 MB - ./generate_garf_training_dataset.py -s intevo -r lu177 -d v3 -n 2e9 - ./generate_garf_training_dataset.py -s nm670 -r tc99m -d v1 ./generate_garf_training_dataset.py -s nm670 -r tc99m -d v2 ./generate_garf_training_dataset.py -s nm670 -r lu177 -d v1 - -## Ideal UI ? - -options -- device: intevo, nm670 -- collimator: lehr, etc -- rad: 177lu, Tc99m etc -- digit + windowing: ?????? -- physics, cuts -- rr=50 ? -- nb particle, root size -- Output: stat + root + simu.json + pth - -Naming: -- intevo_melp_177lu_digit_2_v034.pth -- nm670_3p8_lehr_99tcm_digit_33_v034.pth -- nm670_5p8_lehr_99tcm_digit_33_v034.pth - - generate_training_dataset -s intevo -r 177lu --digit x diff --git a/garf_models/test01_helpers.py b/garf_models/test01_helpers.py index 916b3a1..7e2de0f 100755 --- a/garf_models/test01_helpers.py +++ b/garf_models/test01_helpers.py @@ -6,6 +6,7 @@ from digitizers import * import numpy as np + def init_sim(sim): m = gate.g4_units.m mm = gate.g4_units.mm @@ -22,7 +23,6 @@ def init_sim(sim): return stats - def add_vox_iec(sim, spacing, data_path): mm = gate.g4_units.mm # voxelize_iec_phantom -o data/iec_5mm.mha --spacing 5 @@ -53,42 +53,38 @@ def add_vox_source(sim, rad, activity, data_path): print(f"Total activity is {source.activity / Bq}") return source -def add_source_point_test(sim, rad, planes, activity, angle_tolerance, head_type='arf'): + +def add_source_point_test(sim, rad, planes, activity, angle_tolerance): cm = gate.g4_units.cm source = sim.add_source("GenericSource", "source_point") set_source_rad_energy_spectrum(source, rad) source.particle = "gamma" source.direction.type = "iso" source.position.type = "sphere" - source.position.radius = 1*cm + source.position.radius = 1 * cm source.direction.acceptance_angle.volumes = [p.name for p in planes] source.direction.acceptance_angle.skip_policy = "SkipEvents" source.direction.acceptance_angle.intersection_flag = True source.direction.acceptance_angle.normal_flag = True - source.direction.acceptance_angle.normal_vector = [0, 0, -1] - if head_type == 'spect': - source.direction.acceptance_angle.normal_vector = [1, 0, 0] + source.direction.acceptance_angle.normal_vector = [1, 0, 0] source.direction.acceptance_angle.normal_tolerance = angle_tolerance source.activity = activity - source.attached_to = 'b' + source.attached_to = "b" - b = sim.add_volume('Box', 'b') - #b.translation = [3*cm, 10*cm, 2*cm] - b.translation = [3*cm, 0*cm, 2*cm] - b.material = 'G4_AIR' + b = sim.add_volume("Box", "b") + b.translation = [3 * cm, -10 * cm, 2 * cm] + b.material = "G4_AIR" source = sim.add_source("GenericSource", "source_point2") set_source_rad_energy_spectrum(source, rad) source.particle = "gamma" source.direction.type = "iso" source.position.type = "sphere" - source.position.radius = 1*cm + source.position.radius = 1 * cm source.direction.acceptance_angle.volumes = [p.name for p in planes] source.direction.acceptance_angle.skip_policy = "SkipEvents" source.direction.acceptance_angle.intersection_flag = True source.direction.acceptance_angle.normal_flag = True - source.direction.acceptance_angle.normal_vector = [0, 0, -1] - if head_type == 'spect': - source.direction.acceptance_angle.normal_vector = [1, 0, 0] + source.direction.acceptance_angle.normal_vector = [1, 0, 0] source.direction.acceptance_angle.normal_tolerance = angle_tolerance - source.activity = activity/2 \ No newline at end of file + source.activity = activity / 2 diff --git a/garf_models/test01_main1_intevo_source_point_ref.py b/garf_models/test01_main1_intevo_source_point_ref.py index 3112725..6933dc6 100755 --- a/garf_models/test01_main1_intevo_source_point_ref.py +++ b/garf_models/test01_main1_intevo_source_point_ref.py @@ -2,6 +2,7 @@ # -*- coding: utf-8 -*- import opengate.contrib.spect.siemens_intevo as intevo +from opengate.contrib.spect.spect_helpers import add_fake_table from pathlib import Path from digitizers import * from test01_helpers import init_sim, add_source_point_test @@ -13,7 +14,7 @@ def go(): sim = gate.Simulation() # main options - #sim.visu = True + # sim.visu = True sim.visu_type = "qt" sim.random_seed = "auto" sim.number_of_threads = 4 @@ -47,23 +48,27 @@ def go(): proj1 = add_intevo_digitizer_lu177_v3( sim, crystal1, f"digitizer1", spectrum_channel=False ) - proj1.output_filename = "projection_1.mhd" - head2, colli2, crystal2 = intevo.add_spect_head( sim, "spect2", collimator_type=colli_type, debug=sim.visu == True ) proj2 = add_intevo_digitizer_lu177_v3( sim, crystal2, f"digitizer2", spectrum_channel=False ) + + # output names + proj1.output_filename = "projection_1.mhd" proj2.output_filename = "projection_2.mhd" + # for visualisation debug + # table = add_fake_table(sim) + # table.translation = [0, 320, 0] + # rotate - nb_angle = 3 - intevo.rotate_gantry(head1, radius, 0, initial_rotation='spect') - intevo.rotate_gantry(head2, radius, 180, initial_rotation='spect') + intevo.rotate_gantry(head1, radius, 0) + intevo.rotate_gantry(head2, radius, 180) # source for test - add_source_point_test(sim, rad, (head1, head2), activity, angle_tolerance, head_type='spect') + add_source_point_test(sim, rad, (head1, head2), activity, angle_tolerance) # go sim.run() diff --git a/garf_models/test01_main2_intevo_source_point_arf.py b/garf_models/test01_main2_intevo_source_point_arf.py index bce0b4c..8eebaae 100755 --- a/garf_models/test01_main2_intevo_source_point_arf.py +++ b/garf_models/test01_main2_intevo_source_point_arf.py @@ -13,7 +13,7 @@ def go(): sim = gate.Simulation() # main options - #sim.visu = True + # sim.visu = True sim.visu_type = "qt" sim.random_seed = "auto" sim.number_of_threads = 4 @@ -30,13 +30,13 @@ def go(): radius = 28 * cm rad = "lu177" colli_type = "melp" - activity = 1e8* Bq + activity = 1e8 * Bq angle_tolerance = 10 * deg # visu if sim.visu: sim.number_of_threads = 1 - activity = 1000* Bq + activity = 1000 * Bq # world etc stats = init_sim(sim) @@ -46,19 +46,22 @@ def go(): size = [128 * 2, 128 * 2] pth = Path("pth") / "intevo_lu177_v3.pth" det_plane1, arf1 = intevo.add_arf_detector( - sim, radius, 0, size, spacing, colli_type, "detector", 1, pth + sim, "det1", colli_type, size, spacing, pth ) det_plane2, arf2 = intevo.add_arf_detector( - sim, radius, 180, size, spacing, colli_type, "detector", 2, pth + sim, "det2", colli_type, size, spacing, pth ) + + # output names + arf1.output_filename = "projection_1.mhd" + arf2.output_filename = "projection_2.mhd" + # compute the gantry rotations - nb_angle = 3 - intevo.rotate_gantry(det_plane1, radius, 0, initial_rotation="arf") - intevo.rotate_gantry(det_plane2, radius, 180, initial_rotation="arf") - det_planes = [det_plane1, det_plane2] + intevo.rotate_gantry(det_plane1, radius, 0) + intevo.rotate_gantry(det_plane2, radius, 180) # source for test - add_source_point_test(sim, rad, det_planes, activity, angle_tolerance) + add_source_point_test(sim, rad, (det_plane1, det_plane2), activity, angle_tolerance) # go sim.run() diff --git a/garf_models/test01_main3_intevo_source_point_free_flight.py b/garf_models/test01_main3_intevo_source_point_free_flight.py new file mode 100755 index 0000000..654edfb --- /dev/null +++ b/garf_models/test01_main3_intevo_source_point_free_flight.py @@ -0,0 +1,78 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate.contrib.spect.siemens_intevo as intevo +from pathlib import Path +from digitizers import * +from test01_helpers import init_sim, add_source_point_test + + +def go(): + + # create the simulation + sim = gate.Simulation() + + # main options + # sim.visu = True + sim.visu_type = "qt" + sim.random_seed = "auto" + sim.number_of_threads = 4 + sim.progress_bar = True + sim.output_dir = Path("test01") / f"free_flight" + + # units + mm = gate.g4_units.mm + cm = gate.g4_units.cm + deg = gate.g4_units.deg + Bq = gate.g4_units.Bq + + # options + radius = 28 * cm + rad = "lu177" + colli_type = "melp" + activity = 1e8 * Bq + angle_tolerance = 10 * deg + + # visu + if sim.visu: + sim.number_of_threads = 1 + activity = 1000 * Bq + + # world etc + stats = init_sim(sim) + + # set the two spect heads + spacing = [4.7951998710632 * mm / 2, 4.7951998710632 * mm / 2] + size = [128 * 2, 128 * 2] + pth = Path("pth") / "intevo_lu177_v3.pth" + det_plane1, arf1 = intevo.add_arf_detector( + sim, "det1", colli_type, size, spacing, pth + ) + det_plane2, arf2 = intevo.add_arf_detector( + sim, "det2", colli_type, size, spacing, pth + ) + + # output names + arf1.output_filename = "projection_1.mhd" + arf2.output_filename = "projection_2.mhd" + + # compute the gantry rotations + intevo.rotate_gantry(det_plane1, radius, 0) + intevo.rotate_gantry(det_plane2, radius, 180) + + # source for test + add_source_point_test(sim, rad, (det_plane1, det_plane2), activity, angle_tolerance) + + # FF + ff = sim.add_actor("FreeFlightActor", "ff") + ff.attached_to = "world" + ff.particles = "gamma" + + # go + sim.run() + print(stats) + + +# -------------------------------------------------------------------------- +if __name__ == "__main__": + go() diff --git a/garf_models/test01_main4_compare.py b/garf_models/test01_main4_compare.py index 2d2a5a8..5de8988 100755 --- a/garf_models/test01_main4_compare.py +++ b/garf_models/test01_main4_compare.py @@ -1,16 +1,13 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- -import opengate.contrib.spect.siemens_intevo as intevo from pathlib import Path -from digitizers import * -from test01_helpers import init_sim, add_source_point_test from opengate.tests import utility def go(): - folder_base = Path('test01') + folder_base = Path("test01") folder_ref = folder_base / "reference" folder_arf = folder_base / "arf" folder_ff = folder_base / "free_flight" @@ -20,28 +17,62 @@ def go(): is_ok = True im_name = "projection_1.mhd" - is_ok = utility.assert_images( - folder_ref / im_name, - folder_arf / im_name, - stats_ref, - tolerance=38, - ignore_value_data1=0, - sum_tolerance=8, - axis="x", - ) and is_ok + is_ok = ( + utility.assert_images( + folder_ref / im_name, + folder_arf / im_name, + stats_ref, + tolerance=38, + ignore_value_data1=0, + sum_tolerance=8, + axis="x", + ) + and is_ok + ) im_name = "projection_2.mhd" - is_ok = utility.assert_images( - folder_ref / im_name, - folder_arf / im_name, - stats_ref, - tolerance=38, - ignore_value_data1=0, - sum_tolerance=8, - axis="x", - ) and is_ok - utility.test_ok(is_ok) + is_ok = ( + utility.assert_images( + folder_ref / im_name, + folder_arf / im_name, + stats_ref, + tolerance=38, + ignore_value_data1=0, + sum_tolerance=8, + axis="x", + ) + and is_ok + ) + + im_name = "projection_1.mhd" + is_ok = ( + utility.assert_images( + folder_ref / im_name, + folder_ff / im_name, + stats_ref, + tolerance=38, + ignore_value_data1=0, + sum_tolerance=8, + axis="x", + ) + and is_ok + ) + im_name = "projection_2.mhd" + is_ok = ( + utility.assert_images( + folder_ref / im_name, + folder_ff / im_name, + stats_ref, + tolerance=38, + ignore_value_data1=0, + sum_tolerance=8, + axis="x", + ) + and is_ok + ) + + utility.test_ok(is_ok) # -------------------------------------------------------------------------- From a6232560fc83de49865b1a62acef9052568efb43 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Fri, 20 Dec 2024 17:09:25 +0100 Subject: [PATCH 04/20] remove unused functions --- garf_models/garf_models_helpers.py | 41 ------------------------------ 1 file changed, 41 deletions(-) diff --git a/garf_models/garf_models_helpers.py b/garf_models/garf_models_helpers.py index 95663cf..663e6ed 100755 --- a/garf_models/garf_models_helpers.py +++ b/garf_models/garf_models_helpers.py @@ -74,30 +74,6 @@ def get_collimator_from_rad(spect, rad): return all_collimators[spect][rad] -def add_spect_imaging_device_OLD_TO_REMOVE(sim, spect, rad, crystal_size): - if spect not in all_spects: - fatal(f'Unknown spect system "{spect}", known are: {all_spects} ') - colli_type = get_collimator_from_rad(spect, rad) - mm = gate.g4_units.mm - head, colli, crystal = None, None, None - if spect == "intevo": - head, colli, crystal = intevo.add_spect_head( - sim, "spect", collimator_type=colli_type, debug=sim.visu == True - ) - intevo.set_head_orientation_OLD_TO_REMOVE(head, colli_type, radius=0 * mm) - if spect == "nm670": - head, colli, crystal = nm670.add_spect_head( - sim, - "spect", - collimator_type=colli_type, - debug=sim.visu == True, - crystal_size=crystal_size, - ) - - print("Head position", head.translation) - return head, colli, crystal, colli_type - - def add_spect_imaging_device(sim, spect, rad, crystal_size): if spect not in all_spects: fatal(f'Unknown spect system "{spect}", known are: {all_spects} ') @@ -158,23 +134,6 @@ def add_digitizer(sim, spect, rad, digitizer, crystal): return ew -def add_arf_OLD_TO_REMOVE(sim, spect, head, colli_type, ew, rr): - detector_plane = None - arf = None - if spect not in all_spects: - fatal(f'Unknown spect system "{spect}", known are: {all_spects} ') - if spect == "nm670": - detector_plane, arf = nm670.add_actor_for_arf_training_dataset( - sim, head, colli_type, ew, rr=rr - ) - if spect == "intevo": - detector_plane, arf = intevo.add_actor_for_arf_training_dataset( - sim, colli_type, ew, rr=rr - ) - print("Plane position", detector_plane.translation) - return detector_plane, arf - - def add_arf_training(sim, spect, head, colli_type, ew, rr): detector_plane = None arf = None From bd5560f87173864453d40445eac5144c1c91fa4b Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Sat, 4 Jan 2025 14:50:37 +0100 Subject: [PATCH 05/20] update to correct arf 036 --- ...arf_v034_test.json => train_arf_v035.json} | 2 +- garf_models/pth/train_arf_v036.json | 33 +++++++ garf_models/readme.md | 10 +- garf_models/test01_helpers.py | 4 +- .../test01_main2_intevo_source_point_arf.py | 2 +- ...1_main3_intevo_source_point_free_flight.py | 2 +- garf_models/test01_main4_compare.py | 32 +++++-- garf_models/test02_main1_intevo_iec_ref.py | 79 ++++++++++++++++ garf_models/test02_main2_intevo_iec_arf.py | 87 ++++++++++++++++++ garf_models/test02_main3_intevo_iec_ff.py | 91 +++++++++++++++++++ garf_models/test02_main4_compare.py | 86 ++++++++++++++++++ 11 files changed, 414 insertions(+), 14 deletions(-) rename garf_models/pth/{train_arf_v034_test.json => train_arf_v035.json} (95%) create mode 100644 garf_models/pth/train_arf_v036.json create mode 100755 garf_models/test02_main1_intevo_iec_ref.py create mode 100755 garf_models/test02_main2_intevo_iec_arf.py create mode 100755 garf_models/test02_main3_intevo_iec_ff.py create mode 100755 garf_models/test02_main4_compare.py diff --git a/garf_models/pth/train_arf_v034_test.json b/garf_models/pth/train_arf_v035.json similarity index 95% rename from garf_models/pth/train_arf_v034_test.json rename to garf_models/pth/train_arf_v035.json index 0b975ef..b992088 100644 --- a/garf_models/pth/train_arf_v034_test.json +++ b/garf_models/pth/train_arf_v035.json @@ -12,7 +12,7 @@ "L": 3, "#": "Learning rate", - "learning_rate": 0.00001, + "learning_rate": 0.0001, "#": "Max nb of epoch (iteration)", "epoch_max": 100, diff --git a/garf_models/pth/train_arf_v036.json b/garf_models/pth/train_arf_v036.json new file mode 100644 index 0000000..7706fb0 --- /dev/null +++ b/garf_models/pth/train_arf_v036.json @@ -0,0 +1,33 @@ +{ + "#": "Seed for random generator", + "seed": 123654, + + "#": "Russian Roulette param", + "RR": 50, + + "#": "Number of neurons by layer", + "H": 200, + + "#": "Number of layer", + "L": 3, + + "#": "Learning rate", + "learning_rate": 0.0001, + + "#": "Max nb of epoch (iteration)", + "epoch_max": 200, + "#epoch_max": 4, + + "#": "Nb of epoch without progress to stop earlier", + "early_stopping": 5, + + "#": "Nb of samples by batch", + "batch_size": 2000, + + "#": "Nb of batch per epoch NOT USED !!!!", + "batch_per_epoch": 200, + + "#": "Every epoch: store the current model (in epoch)", + "epoch_store_every": 50000 + +} diff --git a/garf_models/readme.md b/garf_models/readme.md index 7ac449c..cb623d0 100644 --- a/garf_models/readme.md +++ b/garf_models/readme.md @@ -1,13 +1,21 @@ +## test02 + +- idem test01 with voxelized IEC phantom and voxelized source +- ref done 2e4 Bq +- no AA for arf (need scatter) +- AA for ff, direct only + + ## test01 - Intevo - Lu177, melp - no phantom (in Air) - 2 point-sources + AA - +- good one is intevo_lu177_v3_v036.pth ## garf generation diff --git a/garf_models/test01_helpers.py b/garf_models/test01_helpers.py index 7e2de0f..1c4f034 100755 --- a/garf_models/test01_helpers.py +++ b/garf_models/test01_helpers.py @@ -49,8 +49,8 @@ def add_vox_source(sim, rad, activity, data_path): source.particle = "gamma" _, volumes = nemaiec.get_default_sphere_centers_and_volumes() print(f"Volumes are {volumes}") - source.activity = activity * np.array(volumes).sum() - print(f"Total activity is {source.activity / Bq}") + source.activity = activity * np.array(volumes).sum() / sim.number_of_threads + print(f"Total activity is {(source.activity*sim.number_of_threads) / Bq}") return source diff --git a/garf_models/test01_main2_intevo_source_point_arf.py b/garf_models/test01_main2_intevo_source_point_arf.py index 8eebaae..6ac56c7 100755 --- a/garf_models/test01_main2_intevo_source_point_arf.py +++ b/garf_models/test01_main2_intevo_source_point_arf.py @@ -44,7 +44,7 @@ def go(): # set the two spect heads spacing = [4.7951998710632 * mm / 2, 4.7951998710632 * mm / 2] size = [128 * 2, 128 * 2] - pth = Path("pth") / "intevo_lu177_v3.pth" + pth = Path("pth") / "intevo_lu177_v3_v036.pth" det_plane1, arf1 = intevo.add_arf_detector( sim, "det1", colli_type, size, spacing, pth ) diff --git a/garf_models/test01_main3_intevo_source_point_free_flight.py b/garf_models/test01_main3_intevo_source_point_free_flight.py index 654edfb..309a271 100755 --- a/garf_models/test01_main3_intevo_source_point_free_flight.py +++ b/garf_models/test01_main3_intevo_source_point_free_flight.py @@ -44,7 +44,7 @@ def go(): # set the two spect heads spacing = [4.7951998710632 * mm / 2, 4.7951998710632 * mm / 2] size = [128 * 2, 128 * 2] - pth = Path("pth") / "intevo_lu177_v3.pth" + pth = Path("pth") / "intevo_lu177_v3_v036.pth" det_plane1, arf1 = intevo.add_arf_detector( sim, "det1", colli_type, size, spacing, pth ) diff --git a/garf_models/test01_main4_compare.py b/garf_models/test01_main4_compare.py index 5de8988..42a459c 100755 --- a/garf_models/test01_main4_compare.py +++ b/garf_models/test01_main4_compare.py @@ -22,9 +22,9 @@ def go(): folder_ref / im_name, folder_arf / im_name, stats_ref, - tolerance=38, + tolerance=45, ignore_value_data1=0, - sum_tolerance=8, + sum_tolerance=11, axis="x", ) and is_ok @@ -36,23 +36,24 @@ def go(): folder_ref / im_name, folder_arf / im_name, stats_ref, - tolerance=38, + tolerance=45, ignore_value_data1=0, - sum_tolerance=8, + sum_tolerance=9, axis="x", ) and is_ok ) + print() im_name = "projection_1.mhd" is_ok = ( utility.assert_images( folder_ref / im_name, folder_ff / im_name, stats_ref, - tolerance=38, + tolerance=45, ignore_value_data1=0, - sum_tolerance=8, + sum_tolerance=11, axis="x", ) and is_ok @@ -64,9 +65,24 @@ def go(): folder_ref / im_name, folder_ff / im_name, stats_ref, - tolerance=38, + tolerance=45, ignore_value_data1=0, - sum_tolerance=8, + sum_tolerance=9, + axis="x", + ) + and is_ok + ) + + print() + im_name = "projection_1.mhd" + is_ok = ( + utility.assert_images( + folder_arf / im_name, + folder_ff / im_name, + stats_ref, + tolerance=12, + ignore_value_data1=0, + sum_tolerance=1, axis="x", ) and is_ok diff --git a/garf_models/test02_main1_intevo_iec_ref.py b/garf_models/test02_main1_intevo_iec_ref.py new file mode 100755 index 0000000..b0292c6 --- /dev/null +++ b/garf_models/test02_main1_intevo_iec_ref.py @@ -0,0 +1,79 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate.contrib.spect.siemens_intevo as intevo +from pathlib import Path +from digitizers import * +from test01_helpers import add_vox_iec, add_vox_source, init_sim + + +def go(): + + # create the simulation + sim = gate.Simulation() + + # main options + # sim.visu = True + sim.visu_type = "qt" + sim.random_seed = "auto" + sim.number_of_threads = 8 + sim.progress_bar = True + sim.output_dir = Path("test02") / f"reference" + data_path = Path("data") + + # units + cm = gate.g4_units.cm + Bq = gate.g4_units.Bq + + # options + radius = 28 * cm + rad = "lu177" + colli_type = "melp" + activity = 2e4 * Bq + + # visu + if sim.visu: + sim.number_of_threads = 1 + activity = 0.1 * Bq + sim.output_dir = Path("test02") / f"visu" + + # world etc + stats = init_sim(sim) + + # set the spect head + head1, colli1, crystal1 = intevo.add_spect_head( + sim, "spect1", collimator_type=colli_type, debug=sim.visu == True + ) + proj1 = add_intevo_digitizer_lu177_v3( + sim, crystal1, f"digitizer1", spectrum_channel=False + ) + + head2, colli2, crystal2 = intevo.add_spect_head( + sim, "spect2", collimator_type=colli_type, debug=sim.visu == True + ) + proj2 = add_intevo_digitizer_lu177_v3( + sim, crystal2, f"digitizer2", spectrum_channel=False + ) + + # output names + proj1.output_filename = "projection_1.mhd" + proj2.output_filename = "projection_2.mhd" + + # rotate + intevo.rotate_gantry(head1, radius, 0) + intevo.rotate_gantry(head2, radius, 180) + + # add voxelized iec + add_vox_iec(sim, spacing=4, data_path=data_path) + + # add iec voxelized source + add_vox_source(sim, rad, activity, data_path) + + # go + sim.run() + print(stats) + + +# -------------------------------------------------------------------------- +if __name__ == "__main__": + go() diff --git a/garf_models/test02_main2_intevo_iec_arf.py b/garf_models/test02_main2_intevo_iec_arf.py new file mode 100755 index 0000000..d9a9c8f --- /dev/null +++ b/garf_models/test02_main2_intevo_iec_arf.py @@ -0,0 +1,87 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate.contrib.spect.siemens_intevo as intevo +from pathlib import Path +from digitizers import * +from test01_helpers import add_vox_iec, add_vox_source, init_sim + + +def go(): + + # create the simulation + sim = gate.Simulation() + + # main options + # sim.visu = True + sim.visu_type = "qt" + sim.random_seed = "auto" + sim.number_of_threads = 4 + sim.progress_bar = True + sim.output_dir = Path("test02") / f"arf" + data_path = Path("data") + + # units + mm = gate.g4_units.mm + cm = gate.g4_units.cm + m = gate.g4_units.m + Bq = gate.g4_units.Bq + deg = gate.g4_units.deg + + # options + radius = 28 * cm + rad = "lu177" + colli_type = "melp" + activity = 2e3 * Bq + angle_tolerance = 10 * deg + + # visu + if sim.visu: + sim.number_of_threads = 1 + activity = 0.1 * Bq + sim.output_dir = Path("test02") / f"visu" + + # world etc + stats = init_sim(sim) + + # set the two spect heads + spacing = [4.7951998710632 * mm / 2, 4.7951998710632 * mm / 2] + size = [128 * 2, 128 * 2] + pth = Path("pth") / "intevo_lu177_v3_v036.pth" + det_plane1, arf1 = intevo.add_arf_detector( + sim, "det1", colli_type, size, spacing, pth + ) + det_plane2, arf2 = intevo.add_arf_detector( + sim, "det2", colli_type, size, spacing, pth + ) + + # output names + arf1.output_filename = "projection_1.mhd" + arf2.output_filename = "projection_2.mhd" + + # compute the gantry rotations + intevo.rotate_gantry(det_plane1, radius, 0) + intevo.rotate_gantry(det_plane2, radius, 180) + det_planes = (det_plane1, det_plane2) + + # add voxelized iec + add_vox_iec(sim, spacing=4, data_path=data_path) + + # add iec voxelized source + source = add_vox_source(sim, rad, activity, data_path) + """source.direction.acceptance_angle.volumes = [h.name for h in det_planes] + source.direction.acceptance_angle.skip_policy = "SkipEvents" + source.direction.acceptance_angle.intersection_flag = True + source.direction.acceptance_angle.normal_flag = True + source.direction.acceptance_angle.normal_vector = [1, 0, 0] + source.direction.acceptance_angle.normal_tolerance = angle_tolerance + """ + + # go + sim.run() + print(stats) + + +# -------------------------------------------------------------------------- +if __name__ == "__main__": + go() diff --git a/garf_models/test02_main3_intevo_iec_ff.py b/garf_models/test02_main3_intevo_iec_ff.py new file mode 100755 index 0000000..2d3edca --- /dev/null +++ b/garf_models/test02_main3_intevo_iec_ff.py @@ -0,0 +1,91 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import opengate.contrib.spect.siemens_intevo as intevo +from pathlib import Path +from digitizers import * +from test01_helpers import add_vox_iec, add_vox_source, init_sim + + +def go(): + + # create the simulation + sim = gate.Simulation() + + # main options + # sim.visu = True + sim.visu_type = "qt" + sim.random_seed = "auto" + sim.number_of_threads = 4 + sim.progress_bar = True + sim.output_dir = Path("test02") / f"free_flight" + data_path = Path("data") + + # units + mm = gate.g4_units.mm + cm = gate.g4_units.cm + m = gate.g4_units.m + Bq = gate.g4_units.Bq + deg = gate.g4_units.deg + + # options + radius = 28 * cm + rad = "lu177" + colli_type = "melp" + activity = 2e3 * Bq + angle_tolerance = 10 * deg + + # visu + if sim.visu: + sim.number_of_threads = 1 + activity = 0.1 * Bq + sim.output_dir = Path("test02") / f"visu" + + # world etc + stats = init_sim(sim) + + # set the two spect heads + spacing = [4.7951998710632 * mm / 2, 4.7951998710632 * mm / 2] + size = [128 * 2, 128 * 2] + pth = Path("pth") / "intevo_lu177_v3_v036.pth" + det_plane1, arf1 = intevo.add_arf_detector( + sim, "det1", colli_type, size, spacing, pth + ) + det_plane2, arf2 = intevo.add_arf_detector( + sim, "det2", colli_type, size, spacing, pth + ) + + # output names + arf1.output_filename = "projection_1.mhd" + arf2.output_filename = "projection_2.mhd" + + # compute the gantry rotations + intevo.rotate_gantry(det_plane1, radius, 0) + intevo.rotate_gantry(det_plane2, radius, 180) + det_planes = (det_plane1, det_plane2) + + # add voxelized iec + add_vox_iec(sim, spacing=4, data_path=data_path) + + # add iec voxelized source + source = add_vox_source(sim, rad, activity, data_path) + source.direction.acceptance_angle.volumes = [h.name for h in det_planes] + source.direction.acceptance_angle.skip_policy = "SkipEvents" + source.direction.acceptance_angle.intersection_flag = True + source.direction.acceptance_angle.normal_flag = True + source.direction.acceptance_angle.normal_vector = [1, 0, 0] + source.direction.acceptance_angle.normal_tolerance = angle_tolerance + + # FF + ff = sim.add_actor("FreeFlightActor", "ff") + ff.attached_to = "iec" + ff.particles = "gamma" + + # go + sim.run() + print(stats) + + +# -------------------------------------------------------------------------- +if __name__ == "__main__": + go() diff --git a/garf_models/test02_main4_compare.py b/garf_models/test02_main4_compare.py new file mode 100755 index 0000000..3f8c7cd --- /dev/null +++ b/garf_models/test02_main4_compare.py @@ -0,0 +1,86 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +from pathlib import Path +from opengate.tests import utility + + +def go(): + + folder_base = Path("test02") + folder_ref = folder_base / "reference" + folder_arf = folder_base / "arf" + folder_ff = folder_base / "free_flight" + + stats_ref = utility.read_stat_file(folder_ref / "stats.txt") + + is_ok = True + scaling = 10 + + im_name = "projection_1.mhd" + is_ok = ( + utility.assert_images( + folder_ref / im_name, + folder_arf / im_name, + stats_ref, + tolerance=45, + ignore_value_data1=0, + sum_tolerance=11, + scaleImageValuesFactor=scaling, + axis="x", + ) + and is_ok + ) + + im_name = "projection_2.mhd" + is_ok = ( + utility.assert_images( + folder_ref / im_name, + folder_arf / im_name, + stats_ref, + tolerance=45, + ignore_value_data1=0, + sum_tolerance=11, + scaleImageValuesFactor=scaling, + axis="x", + ) + and is_ok + ) + + print() + im_name = "projection_1.mhd" + is_ok = ( + utility.assert_images( + folder_ref / im_name, + folder_ff / im_name, + stats_ref, + tolerance=45, + ignore_value_data1=0, + sum_tolerance=11, + scaleImageValuesFactor=scaling, + axis="x", + ) + and is_ok + ) + + """im_name = "projection_2.mhd" + is_ok = ( + utility.assert_images( + folder_ref / im_name, + folder_ff / im_name, + stats_ref, + tolerance=45, + ignore_value_data1=0, + sum_tolerance=11, + scaleImageValuesFactor=scaling, + axis="x", + ) + and is_ok + )""" + + utility.test_ok(is_ok) + + +# -------------------------------------------------------------------------- +if __name__ == "__main__": + go() From 93cc45fdbf06a98ac4099e331077cfc203d373d6 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Tue, 7 Jan 2025 16:19:10 +0100 Subject: [PATCH 06/20] compare ref and ARF and FF --- garf/bin/garf_compare_image_profile.py | 110 ++++++++++++------------- garf/helpers.py | 60 ++++++++++++++ garf_models/readme.md | 6 +- garf_models/test01_main4_compare.py | 34 ++++---- garf_models/test02_main4_compare.py | 81 ++++++++++++------ 5 files changed, 191 insertions(+), 100 deletions(-) diff --git a/garf/bin/garf_compare_image_profile.py b/garf/bin/garf_compare_image_profile.py index 17accb7..51bc8cd 100755 --- a/garf/bin/garf_compare_image_profile.py +++ b/garf/bin/garf_compare_image_profile.py @@ -14,106 +14,98 @@ @click.argument("image1_mhd") @click.argument("image2_mhd") @click.option( - "--events", - "-e", - default=float(1), + "--scaling", + "-s", + default=1.0, help="Scale the image2 by this value before comparing", ) -@click.option("--islice", "-s", default=int(64), help="Image slice for the profile") +@click.option( + "--islice", "-i", default=None, help="Image slice for the profile (middle if None)" +) @click.option("--wslice", "-w", default=int(3), help="Slice width (to smooth)") -def garf_compare_image_profile(image1_mhd, image2_mhd, islice, events, wslice): +@click.option("--output", "-o", default="output.pdf", help="output") +@click.option( + "-rmf", + default=False, + is_flag=True, + help="Remove first slice of ref data (hit slice)", +) +def garf_compare_image_profile( + image1_mhd, image2_mhd, islice, scaling, wslice, rmf, output +): # Load image img_ref = sitk.ReadImage(image1_mhd) img = sitk.ReadImage(image2_mhd) - events = float(events) - islice = int(islice) + scaling = float(scaling) wslice = int(wslice) - # Scale data to the ref nb of particles - img = img * events + # slice + if islice is None: + islice = int(img.GetSize()[0] / 2) + else: + islice = int(islice) # Get the pixels values as np array data_ref = sitk.GetArrayFromImage(img_ref).astype(float) data = sitk.GetArrayFromImage(img).astype(float) + # Scale data to the ref nb of particles + data = data * scaling + + print(f"Reference image shape : {data_ref.shape}") + print(f"Test image shape : {data.shape}") + # Sometimes not same nb of slices -> crop the data_ref if len(data_ref) > len(data): data_ref = data_ref[0 : len(data), :, :] + # Remove first slice ? + if rmf: + data_ref = data_ref[1:, :, :] + # Criterion1: global counts in every windows s_ref = np.sum(data_ref, axis=(1, 2)) s = np.sum(data, axis=(1, 2)) + ratio = (s - s_ref) / s_ref * 100.0 - print("Ref: Singles/Scatter/Peak1/Peak2: {}".format(s_ref)) - print("Img: WARNING/Scatter/Peak1/Peak2: {}".format(s)) - print( - "% diff : WARNING/Scatter/Peak1/Peak2: {}".format((s - s_ref) / s_ref * 100.0) - ) + # global counts + print(f"Global counts, reference : {s_ref}") + print(f"Global counts, test image: {s}") + print(f"Global counts, % diff : {ratio} %") # Profiles - # data image: !eee!Z,Y,X p_ref = np.mean(data_ref[:, islice - wslice : islice + wslice - 1, :], axis=1) p = np.mean(data[:, islice - wslice : islice + wslice - 1, :], axis=1) - x = np.arange(0, 128, 1) + x = np.arange(0, data.shape[1], 1) + + # max + vmax_ref = np.max(p_ref[1:, :]) + vmax = np.max(p[1:, :]) + print(f"Max value in ref image : {vmax_ref}") + print(f"Max value in test image : {vmax}") + # nb of energy windows nb_ene = len(data) print("Nb of energy windows: ", nb_ene) + win = [f"win {i}" for i in np.arange(nb_ene)] - if nb_ene == 3: # Tc99m - win = ["WARNING", "Scatter", "Peak 140"] - - if nb_ene == 6: # In111 - win = ["WARNING", "Scatter1", "Peak171", "Scatter2", "Scatter3", "Peak245"] - - if nb_ene == 7: # Lu177 - win = [ - "WARNING", - "Scatter1", - "Peak113", - "Scatter2", - "Scatter3", - "Peak208", - "Scatter4", - ] - - if nb_ene == 8: - win = [ - "WARNING", - "Scatter1", - "Peak364", - "Scatter2", - "Scatter3", - "Scatter4", - "Peak637", - "Peak722", - ] - - fig, ax = plt.subplots(ncols=nb_ene - 1, nrows=1, figsize=(35, 5)) - - i = 1 - vmax = np.max(p_ref[1:, :]) - vmax = np.max(p[1:, :]) - print("Max value in ref image for the scale : {}".format(vmax)) - + # figure + fig, ax = plt.subplots(ncols=nb_ene, nrows=1, figsize=(35, 5)) fs = 12 - plt.rc("font", size=fs) - while i < nb_ene: - a = ax[i - 1] - + for i in range(nb_ene): + a = ax[i] a.plot(x, p_ref[i], "g", label="Analog", alpha=0.5, linewidth=2.0) a.plot(x, p[i], "k--", label="ARF", alpha=0.9, linewidth=1.0) a.set_title(win[i], fontsize=fs + 5) a.legend(loc="best") - # a.labelsize = 40 a.tick_params(labelsize=fs) - # a.set_ylim([0, vmax]) i += 1 plt.suptitle("Compare " + image1_mhd + " vs " + image2_mhd + " w=" + str(wslice)) plt.tight_layout() plt.subplots_adjust(top=0.85) - plt.savefig("output.pdf") + plt.savefig(output) plt.show() diff --git a/garf/helpers.py b/garf/helpers.py index 8ef75fa..8a48d44 100644 --- a/garf/helpers.py +++ b/garf/helpers.py @@ -3,6 +3,8 @@ import os import uproot import torch +import SimpleITK as sitk +import matplotlib.pyplot as plt def load_training_dataset(filename): @@ -225,3 +227,61 @@ def get_gpu_device(gpu_mode): return get_gpu_device("cpu") return current_gpu_mode, current_gpu_device + + +def plot_spect_projection( + image1_mhd, image2_mhd, scaling=1, islice=None, wslice=1, win_labels=None +): + # Load image + img_ref = sitk.ReadImage(image1_mhd) + img = sitk.ReadImage(image2_mhd) + scaling = float(scaling) + wslice = int(wslice) + + # slice + if islice is None: + islice = int(img.GetSize()[0] / 2) + else: + islice = int(islice) + + # Get the pixels values as np array + data_ref = sitk.GetArrayFromImage(img_ref).astype(float) + data = sitk.GetArrayFromImage(img).astype(float) + + # Scale data to the ref nb of particles + data = data * scaling + + # Profiles + p_ref = np.mean(data_ref[:, islice - wslice : islice + wslice - 1, :], axis=1) + p = np.mean(data[:, islice - wslice : islice + wslice - 1, :], axis=1) + x = np.arange(0, data.shape[1], 1) + + # nb of energy windows + nb_ene = len(data) + if win_labels is None: + win_labels = [f"win {i}" for i in np.arange(nb_ene)] + + # Criterion1: global counts in every windows + s_ref = np.sum(data_ref, axis=(1, 2)) + s = np.sum(data, axis=(1, 2)) + ratio = (s - s_ref) / s_ref * 100.0 + + # figure + fig, ax = plt.subplots(ncols=nb_ene, nrows=1, figsize=(35, 5)) + fs = 12 + plt.rc("font", size=fs) + for i in range(nb_ene): + a = ax[i] + a.plot(x, p_ref[i], "g", label="Analog", alpha=0.5, linewidth=2.0) + a.plot(x, p[i], "k--", label="ARF", alpha=0.9, linewidth=1.0) + a.set_title(win_labels[i], fontsize=fs + 5) + a.legend(loc="best") + a.tick_params(labelsize=fs) + i += 1 + + plt.suptitle( + f"Slice ={islice}, w={wslice} {image1_mhd} {image2_mhd} global diff = {ratio}" + ) + plt.tight_layout() + plt.subplots_adjust(top=0.85) + return plt diff --git a/garf_models/readme.md b/garf_models/readme.md index cb623d0..9af7357 100644 --- a/garf_models/readme.md +++ b/garf_models/readme.md @@ -5,8 +5,8 @@ - idem test01 with voxelized IEC phantom and voxelized source - ref done 2e4 Bq -- no AA for arf (need scatter) -- AA for ff, direct only +- no AA for arf (need scatter) -> fit ok +- AA for ff, direct only -> fit not ok because of scatter ## test01 @@ -14,7 +14,7 @@ - Intevo - Lu177, melp - no phantom (in Air) -- 2 point-sources + AA +- 2 point sources + AA - good one is intevo_lu177_v3_v036.pth ## garf generation diff --git a/garf_models/test01_main4_compare.py b/garf_models/test01_main4_compare.py index 42a459c..1934a8d 100755 --- a/garf_models/test01_main4_compare.py +++ b/garf_models/test01_main4_compare.py @@ -11,79 +11,83 @@ def go(): folder_ref = folder_base / "reference" folder_arf = folder_base / "arf" folder_ff = folder_base / "free_flight" - stats_ref = utility.read_stat_file(folder_ref / "stats.txt") - is_ok = True + print(f"Ref vs ARF (projection1) =======================") im_name = "projection_1.mhd" is_ok = ( utility.assert_images( folder_ref / im_name, folder_arf / im_name, stats_ref, - tolerance=45, - ignore_value_data1=0, - sum_tolerance=11, + sum_tolerance=9, axis="x", + test_sad=False, + sad_profile_tolerance=10, ) and is_ok ) + print() + print(f"Ref vs ARF (projection2) =======================") im_name = "projection_2.mhd" is_ok = ( utility.assert_images( folder_ref / im_name, folder_arf / im_name, stats_ref, - tolerance=45, - ignore_value_data1=0, + test_sad=False, sum_tolerance=9, axis="x", + sad_profile_tolerance=10, ) and is_ok ) print() + print(f"Ref vs FF (projection1) =======================") im_name = "projection_1.mhd" is_ok = ( utility.assert_images( folder_ref / im_name, folder_ff / im_name, stats_ref, - tolerance=45, - ignore_value_data1=0, - sum_tolerance=11, + test_sad=False, + sum_tolerance=9, axis="x", + sad_profile_tolerance=10, ) and is_ok ) - + print() + print(f"Ref vs FF (projection2) =======================") im_name = "projection_2.mhd" is_ok = ( utility.assert_images( folder_ref / im_name, folder_ff / im_name, stats_ref, - tolerance=45, - ignore_value_data1=0, + test_sad=False, sum_tolerance=9, axis="x", + sad_profile_tolerance=10, ) and is_ok ) print() + print(f"ARF vs FF (projection2) =======================") im_name = "projection_1.mhd" is_ok = ( utility.assert_images( folder_arf / im_name, folder_ff / im_name, stats_ref, - tolerance=12, - ignore_value_data1=0, + test_sad=False, sum_tolerance=1, axis="x", + sad_profile_tolerance=3, ) and is_ok ) diff --git a/garf_models/test02_main4_compare.py b/garf_models/test02_main4_compare.py index 3f8c7cd..b564355 100755 --- a/garf_models/test02_main4_compare.py +++ b/garf_models/test02_main4_compare.py @@ -3,6 +3,7 @@ from pathlib import Path from opengate.tests import utility +from garf.helpers import plot_spect_projection def go(): @@ -11,72 +12,106 @@ def go(): folder_ref = folder_base / "reference" folder_arf = folder_base / "arf" folder_ff = folder_base / "free_flight" - stats_ref = utility.read_stat_file(folder_ref / "stats.txt") is_ok = True scaling = 10 + # ref vs ARF proj1 + print(f"Ref vs ARF (projection1) =======================") im_name = "projection_1.mhd" + f1 = folder_ref / im_name + f2 = folder_arf / im_name is_ok = ( utility.assert_images( - folder_ref / im_name, - folder_arf / im_name, + f1, + f2, stats_ref, - tolerance=45, - ignore_value_data1=0, - sum_tolerance=11, + test_sad=False, + sum_tolerance=15, scaleImageValuesFactor=scaling, axis="x", + sad_profile_tolerance=20, ) and is_ok ) + o = folder_arf / "compare_1.pdf" + plt = plot_spect_projection(f1, f2, scaling=10, islice=117, wslice=3) + plt.savefig(o) + print(f"Profile figure : ", o) + + # ref vs ARF proj2 + print() + print(f"Ref vs ARF (projection2) =======================") im_name = "projection_2.mhd" + f1 = folder_ref / im_name + f2 = folder_arf / im_name is_ok = ( utility.assert_images( - folder_ref / im_name, - folder_arf / im_name, + f1, + f2, stats_ref, - tolerance=45, - ignore_value_data1=0, - sum_tolerance=11, + test_sad=False, + sum_tolerance=15, scaleImageValuesFactor=scaling, axis="x", + sad_profile_tolerance=20, ) and is_ok ) + o = folder_arf / "compare_2.pdf" + plt = plot_spect_projection(f1, f2, scaling=10, islice=117, wslice=3) + plt.savefig(o) + print(f"Profile figure : ", o) + + # ref vs FF proj1 print() + print(f"Ref vs FF (projection1) =======================") im_name = "projection_1.mhd" + f1 = folder_ref / im_name + f2 = folder_ff / im_name is_ok = ( utility.assert_images( - folder_ref / im_name, - folder_ff / im_name, + f1, + f2, stats_ref, - tolerance=45, - ignore_value_data1=0, - sum_tolerance=11, + test_sad=False, + sum_tolerance=70, scaleImageValuesFactor=scaling, axis="x", ) and is_ok ) - """im_name = "projection_2.mhd" + o = folder_ff / "compare_1.pdf" + plt = plot_spect_projection(f1, f2, scaling=10, islice=117, wslice=3) + plt.savefig(o) + print(f"Profile figure : ", o) + + print() + print(f"Ref vs FF (projection1) =======================") + im_name = "projection_2.mhd" + f1 = folder_ref / im_name + f2 = folder_ff / im_name is_ok = ( utility.assert_images( - folder_ref / im_name, - folder_ff / im_name, + f1, + f2, stats_ref, - tolerance=45, - ignore_value_data1=0, - sum_tolerance=11, + test_sad=False, + sum_tolerance=70, scaleImageValuesFactor=scaling, axis="x", ) and is_ok - )""" + ) + + o = folder_ff / "compare_2.pdf" + plt = plot_spect_projection(f1, f2, scaling=10, islice=117, wslice=3) + plt.savefig(o) + print(f"Profile figure : ", o) utility.test_ok(is_ok) From cb98f6d61563711bd98ef2277b767c02f9a077f5 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Fri, 27 Jun 2025 18:16:26 +0200 Subject: [PATCH 07/20] old model --- garf/bin/garf_nn_info.py | 4 +- garf/garf_train.py | 198 +++++++++++++++++- garf_models/digitizers.py | 98 --------- garf_models/garf_models_helpers.py | 170 --------------- garf_models/generate_garf_training_dataset.py | 68 ------ garf_models/old/test01_main4_compare.py | 54 ----- garf_models/pth/train_arf_v035.json | 33 --- garf_models/pth/train_arf_v036.json | 33 --- garf_models/readme.md | 32 --- garf_models/test01_helpers.py | 90 -------- .../test01_main1_intevo_source_point_ref.py | 80 ------- .../test01_main2_intevo_source_point_arf.py | 73 ------- ...1_main3_intevo_source_point_free_flight.py | 78 ------- garf_models/test01_main4_compare.py | 100 --------- garf_models/test02_main1_intevo_iec_ref.py | 79 ------- garf_models/test02_main2_intevo_iec_arf.py | 87 -------- garf_models/test02_main3_intevo_iec_ff.py | 91 -------- garf_models/test02_main4_compare.py | 121 ----------- 18 files changed, 200 insertions(+), 1289 deletions(-) delete mode 100755 garf_models/digitizers.py delete mode 100755 garf_models/garf_models_helpers.py delete mode 100755 garf_models/generate_garf_training_dataset.py delete mode 100755 garf_models/old/test01_main4_compare.py delete mode 100644 garf_models/pth/train_arf_v035.json delete mode 100644 garf_models/pth/train_arf_v036.json delete mode 100644 garf_models/readme.md delete mode 100755 garf_models/test01_helpers.py delete mode 100755 garf_models/test01_main1_intevo_source_point_ref.py delete mode 100755 garf_models/test01_main2_intevo_source_point_arf.py delete mode 100755 garf_models/test01_main3_intevo_source_point_free_flight.py delete mode 100755 garf_models/test01_main4_compare.py delete mode 100755 garf_models/test02_main1_intevo_iec_ref.py delete mode 100755 garf_models/test02_main2_intevo_iec_arf.py delete mode 100755 garf_models/test02_main3_intevo_iec_ff.py delete mode 100755 garf_models/test02_main4_compare.py diff --git a/garf/bin/garf_nn_info.py b/garf/bin/garf_nn_info.py index 18256fb..47eef12 100755 --- a/garf/bin/garf_nn_info.py +++ b/garf/bin/garf_nn_info.py @@ -23,7 +23,9 @@ def garf_nn_info(filename_pth): loss_values = p["loss_values"] x = np.arange(0, len(loss_values)) - plt.plot(x, loss_values) + # plt.plot(x, loss_values) + # remove first and last + plt.plot(x[1:-1], loss_values[1:-1]) plt.tight_layout() plt.show() diff --git a/garf/garf_train.py b/garf/garf_train.py index c1bde60..50b6f5a 100644 --- a/garf/garf_train.py +++ b/garf/garf_train.py @@ -76,7 +76,7 @@ def nn_get_optimiser(model_data, model): return optimizer, scheduler -def train_nn(x_train, y_train, params): +def train_nn_OLD_CORRECT(x_train, y_train, params): """ Train the ARF neural network. @@ -259,3 +259,199 @@ def train_nn(x_train, y_train, params): model_data["best_loss"] = best_loss return nn + + +def train_nn(x_train, y_train, params): + """ + Train the ARF neural network. + + x_train -- x samples (3 dimensions: theta, phi, E) + y_train -- output probabilities vector for N energy windows + params -- dictionary of parameters and options + + params contains: + - n_ene_win + - batch_size + - batch_per_epoch + - epoch_store_every + - H + - L + - epoch_max + - early_stopping + - gpu_mode : auto cpu gpu + """ + + # Initialization + x_train, y_train, model_data, N = nn_prepare_data(x_train, y_train, params) + + # One-hot encoding + print("One-hot encoding") + y_vals, y_train = np.unique(y_train, return_inverse=True) + n_ene_win = len(y_vals) + print("Number of energy windows:", n_ene_win) + model_data["n_ene_win"] = n_ene_win + + # Device type + current_gpu_mode, current_gpu_device = get_gpu_device(params["gpu_mode"]) + model_data["current_gpu_mode"] = current_gpu_mode + print(f"Device GPU type is {current_gpu_mode}") + + # Batch parameters + batch_per_epoch = model_data["batch_per_epoch"] + batch_size = model_data["batch_size"] + epoch_store_every = model_data["epoch_store_every"] + + # DataLoader + print("Data loader batch_size", batch_size) + print("Data loader batch_per_epoch", batch_per_epoch) + train_data2 = np.column_stack((x_train, y_train)) + if current_gpu_mode == "mps": + print("With device mps (gpu), convert data to float32", train_data2.dtype) + train_data2 = train_data2.astype(np.float32) + + train_loader2 = DataLoader( + train_data2, + batch_size=batch_size, + num_workers=2, + # pin_memory=True, + # shuffle=True, # if false ~20% faster, seems identical + shuffle=False, # if false ~20% faster, seems identical + drop_last=True, + ) + + # Create the main NN + H = model_data["H"] + L = model_data["L"] + model = Net_v1(H, L, n_ene_win) + + # Create the optimizer + # optimizer, scheduler = nn_get_optimiser(model_data, model) + learning_rate = model_data["learning_rate"] + # optimizer = torch.optim.RMSprop(model.parameters(), lr=learning_rate) + optimizer = torch.optim.Adam( + model.parameters(), + lr=learning_rate, + # weight_decay=1e-4, + ) + # decreasing learning_rate + scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau( + optimizer, "min", verbose=False, patience=3 + ) + + # Main loop initialization + epoch_max = model_data["epoch_max"] + early_stopping = model_data["early_stopping"] + best_loss = np.Inf + best_epoch = 0 + best_train_loss = np.Inf + loss_values = np.zeros(epoch_max + 1) + + # Print parameters + print_nn_params(model_data) + + # create main structures + nn = dict() + nn["model_data"] = model_data + nn["optim"] = dict() + nn["optim"]["model_state"] = [] + nn["optim"]["data"] = [] + previous_best = 9999 + best_epoch_index = 0 + + # set the model to the device (cpu or gpu = cuda or mps) + model.to(current_gpu_device) + + # Main loop + print("\nStart learning ...") + pbar = tqdm(total=epoch_max + 1, disable=not params["progress_bar"]) + epoch = 0 + stop = False + while (not stop) and (epoch < epoch_max): + # Train pass + model.train() + train_loss = 0.0 + n_samples_processed = 0 + + # Loop on the data batch (batch_per_epoch times) + for batch_idx, data in enumerate(train_loader2): + x = data[:, 0:3] + y = data[:, 3] + X = Tensor(x.to(model.fc1.weight.dtype)).to(current_gpu_device) + Y = Tensor(y).to(current_gpu_device).long() + + # Forward pass + Y_out = model(X) + + # Compute expected loss + # combines log_softmax and nll_loss in a single function + loss = F.cross_entropy(Y_out, Y) + + # Backward pass + loss.backward() + + # Parameter update (gradient descent) + optimizer.step() + optimizer.zero_grad() + batch_size = X.shape[0] # important with variable batch sizes + train_loss += loss.data.item() * batch_size + n_samples_processed += batch_size + + # Stop when batch_per_epoch is reach + if batch_idx == params["batch_per_epoch"]: + break + # end for loop train_loader + + # end of train + train_loss /= n_samples_processed + if train_loss < best_train_loss * (1 - 1e-4): + best_train_loss = train_loss + mean_loss = train_loss + + loss_values[epoch] = mean_loss + if mean_loss < best_loss * (1 - 1e-4): + best_loss = mean_loss + best_epoch = epoch + elif epoch - best_epoch > early_stopping: + tqdm.write( + "{} epochs without improvement, early stop.".format(early_stopping) + ) + stop = True + + # scheduler for learning rate + scheduler.step(mean_loss) + + # FIXME WRONG + # Check if need to print and store this epoch + if best_train_loss < previous_best: + tqdm.write("Epoch {} loss is {:.5f}".format(epoch, best_loss)) + previous_best = best_train_loss + + if ( + (epoch != 0 and epoch % epoch_store_every == 0) + or stop + or epoch >= epoch_max - 1 + ): + optim_data = dict() + print("Store weights", epoch) + optim_data["epoch"] = epoch + optim_data["train_loss"] = train_loss + state = copy.deepcopy(model.state_dict()) + nn["optim"]["model_state"].append(state) + nn["optim"]["data"].append(optim_data) + best_epoch_index = epoch + + # update progress bar + pbar.update(1) + epoch = epoch + 1 + + # end for loop + print("Training done. Best = {:.5f} at epoch {:.0f}".format(best_loss, best_epoch)) + + # prepare data to be saved + model_data["loss_values"] = loss_values + model_data["final_epoch"] = epoch + model_data["best_epoch"] = best_epoch + model_data["best_epoch_index"] = best_epoch_index + model_data["best_loss"] = best_loss + + return nn diff --git a/garf_models/digitizers.py b/garf_models/digitizers.py deleted file mode 100755 index a9e2077..0000000 --- a/garf_models/digitizers.py +++ /dev/null @@ -1,98 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -import opengate as gate -from scipy.spatial.transform import Rotation - - -def add_intevo_digitizer_lu177_v3(sim, crystal_name, name, spectrum_channel=True): - - # hits - hits = sim.add_actor("DigitizerHitsCollectionActor", f"hits_{name}") - hits.attached_to = crystal_name - hits.output_filename = "" # No output - hits.attributes = [ - "PostPosition", - "TotalEnergyDeposit", - "PreStepUniqueVolumeID", - "PostStepUniqueVolumeID", - "GlobalTime", - ] - - # singles - singles = sim.add_actor("DigitizerAdderActor", f"singles_{name}") - singles.attached_to = crystal_name - singles.input_digi_collection = hits.name - # sc.policy = "EnergyWeightedCentroidPosition" - singles.policy = "EnergyWinnerPosition" - singles.output_filename = "" # No output - singles.group_volume = None - - # efficiency actor - eff = sim.add_actor("DigitizerEfficiencyActor", f"singles_{name}_eff") - eff.attached_to = crystal_name - eff.input_digi_collection = singles.name - eff.efficiency = 0.86481 # FIXME probably wrong, to evaluate - eff.efficiency = 1.0 - eff.output_filename = "" # No output - - # energy blur - keV = gate.g4_units.keV - MeV = gate.g4_units.MeV - ene_blur = sim.add_actor("DigitizerBlurringActor", f"singles_{name}_eblur") - ene_blur.output_filename = "" - ene_blur.attached_to = crystal_name - ene_blur.input_digi_collection = eff.name - ene_blur.blur_attribute = "TotalEnergyDeposit" - ene_blur.blur_method = "Linear" - ene_blur.blur_resolution = 0.13 - ene_blur.blur_reference_value = 80 * keV - ene_blur.blur_slope = -0.09 * 1 / MeV - - # spatial blurring - mm = gate.g4_units.mm - spatial_blur = sim.add_actor( - "DigitizerSpatialBlurringActor", f"singles_{name}_sblur" - ) - spatial_blur.output_filename = "" - spatial_blur.attached_to = crystal_name - spatial_blur.input_digi_collection = ene_blur.name - spatial_blur.blur_attribute = "PostPosition" - spatial_blur.blur_fwhm = 3.9 * mm - spatial_blur.keep_in_solid_limits = True - - # energy windows - singles_ene_windows = sim.add_actor( - "DigitizerEnergyWindowsActor", f"{name}_energy_window" - ) - channels = [ - {"name": f"spectrum", "min": 3 * keV, "max": 515 * keV}, - {"name": f"scatter1_{name}", "min": 96 * keV, "max": 104 * keV}, - {"name": f"peak113_{name}", "min": 104.52 * keV, "max": 121.48 * keV}, - {"name": f"scatter2_{name}", "min": 122.48 * keV, "max": 133.12 * keV}, - {"name": f"scatter3_{name}", "min": 176.46 * keV, "max": 191.36 * keV}, - {"name": f"peak208_{name}", "min": 192.4 * keV, "max": 223.6 * keV}, - {"name": f"scatter4_{name}", "min": 224.64 * keV, "max": 243.3 * keV}, - ] - if not spectrum_channel: - channels.pop(0) - singles_ene_windows.attached_to = crystal_name - singles_ene_windows.input_digi_collection = spatial_blur.name - singles_ene_windows.channels = channels - singles_ene_windows.output_filename = "" # No output - - # projection - proj = sim.add_actor("DigitizerProjectionActor", f"{name}_projection") - proj.attached_to = crystal_name - proj.input_digi_collections = [x["name"] for x in channels] - proj.spacing = [4.7951998710632 * mm / 2, 4.7951998710632 * mm / 2] - proj.size = [128 * 2, 128 * 2] - proj.output_filename = "projection.mhd" - proj.origin_as_image_center = True - - # plane orientation - proj.detector_orientation_matrix = Rotation.from_euler( - "yx", (90, 90), degrees=True - ).as_matrix() - - return proj diff --git a/garf_models/garf_models_helpers.py b/garf_models/garf_models_helpers.py deleted file mode 100755 index 663e6ed..0000000 --- a/garf_models/garf_models_helpers.py +++ /dev/null @@ -1,170 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -import opengate.contrib.spect.ge_discovery_nm670 as nm670 -import opengate.contrib.spect.siemens_intevo as intevo -from opengate.sources.base import get_rad_gamma_spectrum -from opengate.exception import fatal -from digitizers import * - -all_spects = ("intevo", "nm670") - -all_collimators = { - "intevo": {"tc99m": "lehr", "lu177": "melp", "in111": "melp", "i131": "he"}, - "nm670": {"tc99m": "lehr", "lu177": "megp", "in111": "megp", "i131": "hegp"}, -} - -all_digitizers = { - "intevo": { - "v1": { - "tc99m": intevo.add_digitizer_tc99m, - "lu177": intevo.add_digitizer_lu177, - }, - "v2": {"tc99m": intevo.add_digitizer_tc99m_v2}, - "v3": {"lu177": add_intevo_digitizer_lu177_v3}, - }, - "nm670": { - "v1": { - "tc99m": nm670.add_digitizer_tc99m, - "lu177": nm670.add_digitizer_lu177, - }, - "v2": {"tc99m": nm670.add_digitizer_tc99m_v2}, - }, -} - - -def init_default_garf_simulation(sim, visu, threads): - # main options - sim.visu = visu - sim.visu_type = "qt" - sim.number_of_threads = threads - sim.progress_bar = True - sim.store_json_archive = True - sim.store_input_files = False - sim.json_archive_filename = "simu.json" - - # units - mm = gate.g4_units.mm - m = gate.g4_units.m - - # world - world = sim.world - world.size = [2 * m, 2 * m, 2 * m] - world.material = "G4_AIR" - - # physics - sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option4" - sim.physics_manager.global_production_cuts.all = 1 * mm - - # add stat actor - stats = sim.add_actor("SimulationStatisticsActor", "stats") - stats.output_filename = "stats.txt" - return stats - - -def get_collimator_from_rad(spect, rad): - spect = spect.lower() - rad = rad.lower() - if spect not in all_collimators: - fatal(f'Unknown spect system "{spect}", known are: {all_collimators.keys()} ') - if rad not in all_collimators[spect]: - fatal( - f'Unknown radionuclide "{rad}", known are: {all_collimators[spect].keys()} ' - ) - return all_collimators[spect][rad] - - -def add_spect_imaging_device(sim, spect, rad, crystal_size): - if spect not in all_spects: - fatal(f'Unknown spect system "{spect}", known are: {all_spects} ') - colli_type = get_collimator_from_rad(spect, rad) - mm = gate.g4_units.mm - head, colli, crystal = None, None, None - if spect == "intevo": - head, colli, crystal = intevo.add_spect_head( - sim, "spect", collimator_type=colli_type, debug=sim.visu == True - ) - # intevo.set_head_orientation_OLD_TO_REMOVE(head, colli_type, radius=0 * mm) - intevo.rotate_gantry(head, radius=0, start_angle_deg=0) - if spect == "nm670": - head, colli, crystal = nm670.add_spect_head( - sim, - "spect", - collimator_type=colli_type, - debug=sim.visu == True, - crystal_size=crystal_size, - ) - - print("Head position", head.translation) - return head, colli, crystal, colli_type - - -def add_digitizer(sim, spect, rad, digitizer, crystal): - if spect not in all_digitizers: - fatal( - f'Unknown digitizer spect system "{spect}", known are: {all_digitizers.keys()} ' - ) - digitizers = all_digitizers[spect] - if digitizer not in digitizers: - fatal(f'Unknown digitizer "{digitizer}", known are: {digitizers.keys()} ') - digitizers = digitizers[digitizer] - if rad not in digitizers: - fatal( - f'Unknown digitizer "{digitizer}" for radionuclide "{rad}", known are: {digitizers.keys()} ' - ) - - # create the digitizer - f = digitizers[rad] - - # digitizer - f(sim, crystal, f"digitizer") - - # special case for spectrum channel (to remove) - proj = sim.get_actor(f"digitizer_projection") - proj.write_to_disk = False - ew = sim.get_actor(f"digitizer_energy_window") - channels = ew.channels - c = [] - for channel in channels: - if channel.name != "spectrum": - c.append(channel) - ew.channels = c - proj.input_digi_collections = [c["name"] for c in ew.channels] - - return ew - - -def add_arf_training(sim, spect, head, colli_type, ew, rr): - detector_plane = None - arf = None - if spect not in all_spects: - fatal(f'Unknown spect system "{spect}", known are: {all_spects} ') - if spect == "nm670": - detector_plane, arf = nm670.add_actor_for_arf_training_dataset( - sim, head, colli_type, ew, rr=rr - ) - if spect == "intevo": - detector_plane, arf = intevo.add_actor_for_arf_training_dataset( - sim, colli_type, ew, rr=rr - ) - print("Plane position", detector_plane.translation) - return detector_plane, arf - - -def add_source(sim, spect, rad, activity, detector_plane): - if spect not in all_spects: - fatal(f'Unknown spect system "{spect}", known are: {all_spects} ') - MeV = gate.g4_units.MeV - max_energy = get_rad_gamma_spectrum(rad).energies[-1] * 1.05 - print(f"Max energy for {rad} is {max_energy} MeV") - source = None - if spect == "nm670": - source = nm670.add_source_for_arf_training_dataset( - sim, "src", activity, detector_plane, 0.01 * MeV, max_energy - ) - if spect == "intevo": - source = intevo.add_source_for_arf_training_dataset( - sim, "src", activity, detector_plane, 0.01 * MeV, max_energy - ) - print("Source position", source.position.translation) - return source diff --git a/garf_models/generate_garf_training_dataset.py b/garf_models/generate_garf_training_dataset.py deleted file mode 100755 index fbc6d1d..0000000 --- a/garf_models/generate_garf_training_dataset.py +++ /dev/null @@ -1,68 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -from pathlib import Path -import click -from garf_models_helpers import * - -CONTEXT_SETTINGS = dict(help_option_names=["-h", "--help"]) - - -@click.command(context_settings=CONTEXT_SETTINGS) -@click.option("--spect", "-s", default="intevo", help="SPECT system: intevo, nm670") -@click.option("--rad", "-r", required=True, help="Radionuclide 177lu, 99tcm etc") -@click.option("--digitizer", "-d", required=True, help="Digitizer version : v1 v2 etc") -@click.option("--n", "-n", default=5e8, help="Number of particles") -@click.option("--rr", default=50, help="Russian Roulette for GARF") -@click.option("--visu", is_flag=True, default=False, help="visu only") -@click.option("--threads", "-t", default=4, help="number of threads") -@click.option( - "--crystal_size", "-c", default="5/8", help="For nm670 crystal size 3/8 or 5/8" -) -def go(spect, rad, digitizer, n, rr, visu, threads, crystal_size): - - # set the default simulation param - spect_type = spect - sim = gate.Simulation() - arf_basename = f"{spect_type}_{rad}_{digitizer}" - stats = init_default_garf_simulation(sim, visu, threads) - sim.output_dir = Path("output") / arf_basename - - # set the SPECT system - print(f"SPECT type is : {spect_type}") - head, colli, crystal, colli_type = add_spect_imaging_device( - sim, spect_type, rad, crystal_size - ) - print(f"collimator type is : {colli_type}") - - # set the digitizer and the arf - ew = add_digitizer(sim, spect_type, rad, digitizer, crystal) - - # set the arf plane for training - detector_plane, arf = add_arf_training(sim, spect_type, head, colli_type, ew, rr) - - # visu option - if visu: - n = 1e2 - sim.number_of_threads = 1 - sim.output_dir = Path("output") / f"visu" - - # set the source - Bq = gate.g4_units.Bq - activity = n * Bq / sim.number_of_threads - add_source(sim, spect_type, rad, activity, detector_plane) - - # go - sim.run() - - # print results at the end - print(stats) - - # print cmd line to train GARF - pth_json = Path("pth") / "train_arf_v034.json" - pth = Path("pth") / f"{arf_basename}.pth" - print(f"garf_train {pth_json} {arf.get_output_path()} {pth}") - - -# -------------------------------------------------------------------------- -if __name__ == "__main__": - go() diff --git a/garf_models/old/test01_main4_compare.py b/garf_models/old/test01_main4_compare.py deleted file mode 100755 index 65a997d..0000000 --- a/garf_models/old/test01_main4_compare.py +++ /dev/null @@ -1,54 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -import opengate.contrib.spect.siemens_intevo as intevo -from pathlib import Path -from digitizers import * -from test01_helpers import init_sim, add_source_point_test -from opengate.tests import utility - - -def go(): - - folder_base = Path("../test01") - folder_ref = folder_base / "reference" - folder_arf = folder_base / "arf" - folder_ff = folder_base / "free_flight" - - stats_ref = utility.read_stat_file(folder_ref / "stats.txt") - - is_ok = True - - im_name = "projection_1.mhd" - is_ok = ( - utility.assert_images( - folder_ref / im_name, - folder_arf / im_name, - stats_ref, - tolerance=38, - ignore_value_data1=0, - sum_tolerance=8, - axis="x", - ) - and is_ok - ) - - im_name = "projection_2.mhd" - is_ok = ( - utility.assert_images( - folder_ref / im_name, - folder_arf / im_name, - stats_ref, - tolerance=38, - ignore_value_data1=0, - sum_tolerance=8, - axis="x", - ) - and is_ok - ) - utility.test_ok(is_ok) - - -# -------------------------------------------------------------------------- -if __name__ == "__main__": - go() diff --git a/garf_models/pth/train_arf_v035.json b/garf_models/pth/train_arf_v035.json deleted file mode 100644 index b992088..0000000 --- a/garf_models/pth/train_arf_v035.json +++ /dev/null @@ -1,33 +0,0 @@ -{ - "#": "Seed for random generator", - "seed": 123654, - - "#": "Russian Roulette param", - "RR": 50, - - "#": "Number of neurons by layer", - "H": 200, - - "#": "Number of layer", - "L": 3, - - "#": "Learning rate", - "learning_rate": 0.0001, - - "#": "Max nb of epoch (iteration)", - "epoch_max": 100, - "#epoch_max": 4, - - "#": "Nb of epoch without progress to stop earlier", - "early_stopping": 5, - - "#": "Nb of samples by batch", - "batch_size": 2000, - - "#": "Nb of batch per epoch NOT USED !!!!", - "batch_per_epoch": 200, - - "#": "Every epoch: store the current model (in epoch)", - "epoch_store_every": 50000 - -} diff --git a/garf_models/pth/train_arf_v036.json b/garf_models/pth/train_arf_v036.json deleted file mode 100644 index 7706fb0..0000000 --- a/garf_models/pth/train_arf_v036.json +++ /dev/null @@ -1,33 +0,0 @@ -{ - "#": "Seed for random generator", - "seed": 123654, - - "#": "Russian Roulette param", - "RR": 50, - - "#": "Number of neurons by layer", - "H": 200, - - "#": "Number of layer", - "L": 3, - - "#": "Learning rate", - "learning_rate": 0.0001, - - "#": "Max nb of epoch (iteration)", - "epoch_max": 200, - "#epoch_max": 4, - - "#": "Nb of epoch without progress to stop earlier", - "early_stopping": 5, - - "#": "Nb of samples by batch", - "batch_size": 2000, - - "#": "Nb of batch per epoch NOT USED !!!!", - "batch_per_epoch": 200, - - "#": "Every epoch: store the current model (in epoch)", - "epoch_store_every": 50000 - -} diff --git a/garf_models/readme.md b/garf_models/readme.md deleted file mode 100644 index 9af7357..0000000 --- a/garf_models/readme.md +++ /dev/null @@ -1,32 +0,0 @@ - - - -## test02 - -- idem test01 with voxelized IEC phantom and voxelized source -- ref done 2e4 Bq -- no AA for arf (need scatter) -> fit ok -- AA for ff, direct only -> fit not ok because of scatter - - -## test01 - -- Intevo -- Lu177, melp -- no phantom (in Air) -- 2 point sources + AA -- good one is intevo_lu177_v3_v036.pth - -## garf generation - - # osx = 22 min, 346 MB - ./generate_garf_training_dataset.py -s intevo -r lu177 -d v3 -n 2e9 -t 4 - - ./generate_garf_training_dataset.py -s intevo -r tc99m -d v1 - ./generate_garf_training_dataset.py -s intevo -r lu177 -d v1 - ./generate_garf_training_dataset.py -s intevo -r tc99m -d v2 - - ./generate_garf_training_dataset.py -s nm670 -r tc99m -d v1 - ./generate_garf_training_dataset.py -s nm670 -r tc99m -d v2 - ./generate_garf_training_dataset.py -s nm670 -r lu177 -d v1 - diff --git a/garf_models/test01_helpers.py b/garf_models/test01_helpers.py deleted file mode 100755 index 1c4f034..0000000 --- a/garf_models/test01_helpers.py +++ /dev/null @@ -1,90 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -import opengate.contrib.phantoms.nemaiec as nemaiec -from opengate.sources.base import set_source_rad_energy_spectrum -from digitizers import * -import numpy as np - - -def init_sim(sim): - m = gate.g4_units.m - mm = gate.g4_units.mm - # world - world = sim.world - world.size = [2 * m, 2 * m, 2 * m] - world.material = "G4_AIR" - # physics - sim.physics_manager.physics_list_name = "G4EmStandardPhysics_option3" - sim.physics_manager.global_production_cuts.all = 1 * mm - # add stat actor - stats = sim.add_actor("SimulationStatisticsActor", "stats") - stats.output_filename = f"stats.txt" - return stats - - -def add_vox_iec(sim, spacing, data_path): - mm = gate.g4_units.mm - # voxelize_iec_phantom -o data/iec_5mm.mha --spacing 5 - print(f"Phantom: IEC voxelized {spacing}mm") - mhd_filename = data_path / f"iec_{spacing}mm.mha" - labels_filename = data_path / f"iec_{spacing}mm.json" - iec = sim.add_volume("Image", "iec") - iec.image = mhd_filename - nemaiec.create_material(sim) - iec.set_materials_from_voxelisation(labels_filename) - sim.physics_manager.set_production_cut(iec.name, "all", 2 * mm) - return iec - - -def add_vox_source(sim, rad, activity, data_path): - mm = gate.g4_units.mm - Bq = gate.g4_units.Bq - # voxelize_iec_phantom -o data/iec_1mm.mhd --spacing 1 --output_source data/iec_1mm_activity.mhd -a 1 2 3 4 5 6 - iec_source_filename = data_path / "iec_1mm_activity.mhd" - source = sim.add_source("VoxelSource", "src") - source.image = iec_source_filename - source.position.translation = [0, 35 * mm, 0] - set_source_rad_energy_spectrum(source, rad) - source.particle = "gamma" - _, volumes = nemaiec.get_default_sphere_centers_and_volumes() - print(f"Volumes are {volumes}") - source.activity = activity * np.array(volumes).sum() / sim.number_of_threads - print(f"Total activity is {(source.activity*sim.number_of_threads) / Bq}") - return source - - -def add_source_point_test(sim, rad, planes, activity, angle_tolerance): - cm = gate.g4_units.cm - source = sim.add_source("GenericSource", "source_point") - set_source_rad_energy_spectrum(source, rad) - source.particle = "gamma" - source.direction.type = "iso" - source.position.type = "sphere" - source.position.radius = 1 * cm - source.direction.acceptance_angle.volumes = [p.name for p in planes] - source.direction.acceptance_angle.skip_policy = "SkipEvents" - source.direction.acceptance_angle.intersection_flag = True - source.direction.acceptance_angle.normal_flag = True - source.direction.acceptance_angle.normal_vector = [1, 0, 0] - source.direction.acceptance_angle.normal_tolerance = angle_tolerance - source.activity = activity - source.attached_to = "b" - - b = sim.add_volume("Box", "b") - b.translation = [3 * cm, -10 * cm, 2 * cm] - b.material = "G4_AIR" - - source = sim.add_source("GenericSource", "source_point2") - set_source_rad_energy_spectrum(source, rad) - source.particle = "gamma" - source.direction.type = "iso" - source.position.type = "sphere" - source.position.radius = 1 * cm - source.direction.acceptance_angle.volumes = [p.name for p in planes] - source.direction.acceptance_angle.skip_policy = "SkipEvents" - source.direction.acceptance_angle.intersection_flag = True - source.direction.acceptance_angle.normal_flag = True - source.direction.acceptance_angle.normal_vector = [1, 0, 0] - source.direction.acceptance_angle.normal_tolerance = angle_tolerance - source.activity = activity / 2 diff --git a/garf_models/test01_main1_intevo_source_point_ref.py b/garf_models/test01_main1_intevo_source_point_ref.py deleted file mode 100755 index 6933dc6..0000000 --- a/garf_models/test01_main1_intevo_source_point_ref.py +++ /dev/null @@ -1,80 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -import opengate.contrib.spect.siemens_intevo as intevo -from opengate.contrib.spect.spect_helpers import add_fake_table -from pathlib import Path -from digitizers import * -from test01_helpers import init_sim, add_source_point_test - - -def go(): - - # create the simulation - sim = gate.Simulation() - - # main options - # sim.visu = True - sim.visu_type = "qt" - sim.random_seed = "auto" - sim.number_of_threads = 4 - sim.progress_bar = True - sim.output_dir = Path("test01") / f"reference" - - # units - cm = gate.g4_units.cm - deg = gate.g4_units.deg - Bq = gate.g4_units.Bq - - # options - radius = 28 * cm - rad = "lu177" - colli_type = "melp" - activity = 1e8 * Bq - angle_tolerance = 10 * deg - - # visu - if sim.visu: - sim.number_of_threads = 1 - activity = 1000 * Bq - - # world etc - stats = init_sim(sim) - - # set the spect head - head1, colli1, crystal1 = intevo.add_spect_head( - sim, "spect1", collimator_type=colli_type, debug=sim.visu == True - ) - proj1 = add_intevo_digitizer_lu177_v3( - sim, crystal1, f"digitizer1", spectrum_channel=False - ) - head2, colli2, crystal2 = intevo.add_spect_head( - sim, "spect2", collimator_type=colli_type, debug=sim.visu == True - ) - proj2 = add_intevo_digitizer_lu177_v3( - sim, crystal2, f"digitizer2", spectrum_channel=False - ) - - # output names - proj1.output_filename = "projection_1.mhd" - proj2.output_filename = "projection_2.mhd" - - # for visualisation debug - # table = add_fake_table(sim) - # table.translation = [0, 320, 0] - - # rotate - intevo.rotate_gantry(head1, radius, 0) - intevo.rotate_gantry(head2, radius, 180) - - # source for test - add_source_point_test(sim, rad, (head1, head2), activity, angle_tolerance) - - # go - sim.run() - print(stats) - - -# -------------------------------------------------------------------------- -if __name__ == "__main__": - go() diff --git a/garf_models/test01_main2_intevo_source_point_arf.py b/garf_models/test01_main2_intevo_source_point_arf.py deleted file mode 100755 index 6ac56c7..0000000 --- a/garf_models/test01_main2_intevo_source_point_arf.py +++ /dev/null @@ -1,73 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -import opengate.contrib.spect.siemens_intevo as intevo -from pathlib import Path -from digitizers import * -from test01_helpers import init_sim, add_source_point_test - - -def go(): - - # create the simulation - sim = gate.Simulation() - - # main options - # sim.visu = True - sim.visu_type = "qt" - sim.random_seed = "auto" - sim.number_of_threads = 4 - sim.progress_bar = True - sim.output_dir = Path("test01") / f"arf" - - # units - mm = gate.g4_units.mm - cm = gate.g4_units.cm - deg = gate.g4_units.deg - Bq = gate.g4_units.Bq - - # options - radius = 28 * cm - rad = "lu177" - colli_type = "melp" - activity = 1e8 * Bq - angle_tolerance = 10 * deg - - # visu - if sim.visu: - sim.number_of_threads = 1 - activity = 1000 * Bq - - # world etc - stats = init_sim(sim) - - # set the two spect heads - spacing = [4.7951998710632 * mm / 2, 4.7951998710632 * mm / 2] - size = [128 * 2, 128 * 2] - pth = Path("pth") / "intevo_lu177_v3_v036.pth" - det_plane1, arf1 = intevo.add_arf_detector( - sim, "det1", colli_type, size, spacing, pth - ) - det_plane2, arf2 = intevo.add_arf_detector( - sim, "det2", colli_type, size, spacing, pth - ) - - # output names - arf1.output_filename = "projection_1.mhd" - arf2.output_filename = "projection_2.mhd" - - # compute the gantry rotations - intevo.rotate_gantry(det_plane1, radius, 0) - intevo.rotate_gantry(det_plane2, radius, 180) - - # source for test - add_source_point_test(sim, rad, (det_plane1, det_plane2), activity, angle_tolerance) - - # go - sim.run() - print(stats) - - -# -------------------------------------------------------------------------- -if __name__ == "__main__": - go() diff --git a/garf_models/test01_main3_intevo_source_point_free_flight.py b/garf_models/test01_main3_intevo_source_point_free_flight.py deleted file mode 100755 index 309a271..0000000 --- a/garf_models/test01_main3_intevo_source_point_free_flight.py +++ /dev/null @@ -1,78 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -import opengate.contrib.spect.siemens_intevo as intevo -from pathlib import Path -from digitizers import * -from test01_helpers import init_sim, add_source_point_test - - -def go(): - - # create the simulation - sim = gate.Simulation() - - # main options - # sim.visu = True - sim.visu_type = "qt" - sim.random_seed = "auto" - sim.number_of_threads = 4 - sim.progress_bar = True - sim.output_dir = Path("test01") / f"free_flight" - - # units - mm = gate.g4_units.mm - cm = gate.g4_units.cm - deg = gate.g4_units.deg - Bq = gate.g4_units.Bq - - # options - radius = 28 * cm - rad = "lu177" - colli_type = "melp" - activity = 1e8 * Bq - angle_tolerance = 10 * deg - - # visu - if sim.visu: - sim.number_of_threads = 1 - activity = 1000 * Bq - - # world etc - stats = init_sim(sim) - - # set the two spect heads - spacing = [4.7951998710632 * mm / 2, 4.7951998710632 * mm / 2] - size = [128 * 2, 128 * 2] - pth = Path("pth") / "intevo_lu177_v3_v036.pth" - det_plane1, arf1 = intevo.add_arf_detector( - sim, "det1", colli_type, size, spacing, pth - ) - det_plane2, arf2 = intevo.add_arf_detector( - sim, "det2", colli_type, size, spacing, pth - ) - - # output names - arf1.output_filename = "projection_1.mhd" - arf2.output_filename = "projection_2.mhd" - - # compute the gantry rotations - intevo.rotate_gantry(det_plane1, radius, 0) - intevo.rotate_gantry(det_plane2, radius, 180) - - # source for test - add_source_point_test(sim, rad, (det_plane1, det_plane2), activity, angle_tolerance) - - # FF - ff = sim.add_actor("FreeFlightActor", "ff") - ff.attached_to = "world" - ff.particles = "gamma" - - # go - sim.run() - print(stats) - - -# -------------------------------------------------------------------------- -if __name__ == "__main__": - go() diff --git a/garf_models/test01_main4_compare.py b/garf_models/test01_main4_compare.py deleted file mode 100755 index 1934a8d..0000000 --- a/garf_models/test01_main4_compare.py +++ /dev/null @@ -1,100 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -from pathlib import Path -from opengate.tests import utility - - -def go(): - - folder_base = Path("test01") - folder_ref = folder_base / "reference" - folder_arf = folder_base / "arf" - folder_ff = folder_base / "free_flight" - stats_ref = utility.read_stat_file(folder_ref / "stats.txt") - is_ok = True - - print(f"Ref vs ARF (projection1) =======================") - im_name = "projection_1.mhd" - is_ok = ( - utility.assert_images( - folder_ref / im_name, - folder_arf / im_name, - stats_ref, - sum_tolerance=9, - axis="x", - test_sad=False, - sad_profile_tolerance=10, - ) - and is_ok - ) - - print() - print(f"Ref vs ARF (projection2) =======================") - im_name = "projection_2.mhd" - is_ok = ( - utility.assert_images( - folder_ref / im_name, - folder_arf / im_name, - stats_ref, - test_sad=False, - sum_tolerance=9, - axis="x", - sad_profile_tolerance=10, - ) - and is_ok - ) - - print() - print(f"Ref vs FF (projection1) =======================") - im_name = "projection_1.mhd" - is_ok = ( - utility.assert_images( - folder_ref / im_name, - folder_ff / im_name, - stats_ref, - test_sad=False, - sum_tolerance=9, - axis="x", - sad_profile_tolerance=10, - ) - and is_ok - ) - print() - print(f"Ref vs FF (projection2) =======================") - im_name = "projection_2.mhd" - is_ok = ( - utility.assert_images( - folder_ref / im_name, - folder_ff / im_name, - stats_ref, - test_sad=False, - sum_tolerance=9, - axis="x", - sad_profile_tolerance=10, - ) - and is_ok - ) - - print() - print(f"ARF vs FF (projection2) =======================") - im_name = "projection_1.mhd" - is_ok = ( - utility.assert_images( - folder_arf / im_name, - folder_ff / im_name, - stats_ref, - test_sad=False, - sum_tolerance=1, - axis="x", - sad_profile_tolerance=3, - ) - and is_ok - ) - - utility.test_ok(is_ok) - - -# -------------------------------------------------------------------------- -if __name__ == "__main__": - go() diff --git a/garf_models/test02_main1_intevo_iec_ref.py b/garf_models/test02_main1_intevo_iec_ref.py deleted file mode 100755 index b0292c6..0000000 --- a/garf_models/test02_main1_intevo_iec_ref.py +++ /dev/null @@ -1,79 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -import opengate.contrib.spect.siemens_intevo as intevo -from pathlib import Path -from digitizers import * -from test01_helpers import add_vox_iec, add_vox_source, init_sim - - -def go(): - - # create the simulation - sim = gate.Simulation() - - # main options - # sim.visu = True - sim.visu_type = "qt" - sim.random_seed = "auto" - sim.number_of_threads = 8 - sim.progress_bar = True - sim.output_dir = Path("test02") / f"reference" - data_path = Path("data") - - # units - cm = gate.g4_units.cm - Bq = gate.g4_units.Bq - - # options - radius = 28 * cm - rad = "lu177" - colli_type = "melp" - activity = 2e4 * Bq - - # visu - if sim.visu: - sim.number_of_threads = 1 - activity = 0.1 * Bq - sim.output_dir = Path("test02") / f"visu" - - # world etc - stats = init_sim(sim) - - # set the spect head - head1, colli1, crystal1 = intevo.add_spect_head( - sim, "spect1", collimator_type=colli_type, debug=sim.visu == True - ) - proj1 = add_intevo_digitizer_lu177_v3( - sim, crystal1, f"digitizer1", spectrum_channel=False - ) - - head2, colli2, crystal2 = intevo.add_spect_head( - sim, "spect2", collimator_type=colli_type, debug=sim.visu == True - ) - proj2 = add_intevo_digitizer_lu177_v3( - sim, crystal2, f"digitizer2", spectrum_channel=False - ) - - # output names - proj1.output_filename = "projection_1.mhd" - proj2.output_filename = "projection_2.mhd" - - # rotate - intevo.rotate_gantry(head1, radius, 0) - intevo.rotate_gantry(head2, radius, 180) - - # add voxelized iec - add_vox_iec(sim, spacing=4, data_path=data_path) - - # add iec voxelized source - add_vox_source(sim, rad, activity, data_path) - - # go - sim.run() - print(stats) - - -# -------------------------------------------------------------------------- -if __name__ == "__main__": - go() diff --git a/garf_models/test02_main2_intevo_iec_arf.py b/garf_models/test02_main2_intevo_iec_arf.py deleted file mode 100755 index d9a9c8f..0000000 --- a/garf_models/test02_main2_intevo_iec_arf.py +++ /dev/null @@ -1,87 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -import opengate.contrib.spect.siemens_intevo as intevo -from pathlib import Path -from digitizers import * -from test01_helpers import add_vox_iec, add_vox_source, init_sim - - -def go(): - - # create the simulation - sim = gate.Simulation() - - # main options - # sim.visu = True - sim.visu_type = "qt" - sim.random_seed = "auto" - sim.number_of_threads = 4 - sim.progress_bar = True - sim.output_dir = Path("test02") / f"arf" - data_path = Path("data") - - # units - mm = gate.g4_units.mm - cm = gate.g4_units.cm - m = gate.g4_units.m - Bq = gate.g4_units.Bq - deg = gate.g4_units.deg - - # options - radius = 28 * cm - rad = "lu177" - colli_type = "melp" - activity = 2e3 * Bq - angle_tolerance = 10 * deg - - # visu - if sim.visu: - sim.number_of_threads = 1 - activity = 0.1 * Bq - sim.output_dir = Path("test02") / f"visu" - - # world etc - stats = init_sim(sim) - - # set the two spect heads - spacing = [4.7951998710632 * mm / 2, 4.7951998710632 * mm / 2] - size = [128 * 2, 128 * 2] - pth = Path("pth") / "intevo_lu177_v3_v036.pth" - det_plane1, arf1 = intevo.add_arf_detector( - sim, "det1", colli_type, size, spacing, pth - ) - det_plane2, arf2 = intevo.add_arf_detector( - sim, "det2", colli_type, size, spacing, pth - ) - - # output names - arf1.output_filename = "projection_1.mhd" - arf2.output_filename = "projection_2.mhd" - - # compute the gantry rotations - intevo.rotate_gantry(det_plane1, radius, 0) - intevo.rotate_gantry(det_plane2, radius, 180) - det_planes = (det_plane1, det_plane2) - - # add voxelized iec - add_vox_iec(sim, spacing=4, data_path=data_path) - - # add iec voxelized source - source = add_vox_source(sim, rad, activity, data_path) - """source.direction.acceptance_angle.volumes = [h.name for h in det_planes] - source.direction.acceptance_angle.skip_policy = "SkipEvents" - source.direction.acceptance_angle.intersection_flag = True - source.direction.acceptance_angle.normal_flag = True - source.direction.acceptance_angle.normal_vector = [1, 0, 0] - source.direction.acceptance_angle.normal_tolerance = angle_tolerance - """ - - # go - sim.run() - print(stats) - - -# -------------------------------------------------------------------------- -if __name__ == "__main__": - go() diff --git a/garf_models/test02_main3_intevo_iec_ff.py b/garf_models/test02_main3_intevo_iec_ff.py deleted file mode 100755 index 2d3edca..0000000 --- a/garf_models/test02_main3_intevo_iec_ff.py +++ /dev/null @@ -1,91 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -import opengate.contrib.spect.siemens_intevo as intevo -from pathlib import Path -from digitizers import * -from test01_helpers import add_vox_iec, add_vox_source, init_sim - - -def go(): - - # create the simulation - sim = gate.Simulation() - - # main options - # sim.visu = True - sim.visu_type = "qt" - sim.random_seed = "auto" - sim.number_of_threads = 4 - sim.progress_bar = True - sim.output_dir = Path("test02") / f"free_flight" - data_path = Path("data") - - # units - mm = gate.g4_units.mm - cm = gate.g4_units.cm - m = gate.g4_units.m - Bq = gate.g4_units.Bq - deg = gate.g4_units.deg - - # options - radius = 28 * cm - rad = "lu177" - colli_type = "melp" - activity = 2e3 * Bq - angle_tolerance = 10 * deg - - # visu - if sim.visu: - sim.number_of_threads = 1 - activity = 0.1 * Bq - sim.output_dir = Path("test02") / f"visu" - - # world etc - stats = init_sim(sim) - - # set the two spect heads - spacing = [4.7951998710632 * mm / 2, 4.7951998710632 * mm / 2] - size = [128 * 2, 128 * 2] - pth = Path("pth") / "intevo_lu177_v3_v036.pth" - det_plane1, arf1 = intevo.add_arf_detector( - sim, "det1", colli_type, size, spacing, pth - ) - det_plane2, arf2 = intevo.add_arf_detector( - sim, "det2", colli_type, size, spacing, pth - ) - - # output names - arf1.output_filename = "projection_1.mhd" - arf2.output_filename = "projection_2.mhd" - - # compute the gantry rotations - intevo.rotate_gantry(det_plane1, radius, 0) - intevo.rotate_gantry(det_plane2, radius, 180) - det_planes = (det_plane1, det_plane2) - - # add voxelized iec - add_vox_iec(sim, spacing=4, data_path=data_path) - - # add iec voxelized source - source = add_vox_source(sim, rad, activity, data_path) - source.direction.acceptance_angle.volumes = [h.name for h in det_planes] - source.direction.acceptance_angle.skip_policy = "SkipEvents" - source.direction.acceptance_angle.intersection_flag = True - source.direction.acceptance_angle.normal_flag = True - source.direction.acceptance_angle.normal_vector = [1, 0, 0] - source.direction.acceptance_angle.normal_tolerance = angle_tolerance - - # FF - ff = sim.add_actor("FreeFlightActor", "ff") - ff.attached_to = "iec" - ff.particles = "gamma" - - # go - sim.run() - print(stats) - - -# -------------------------------------------------------------------------- -if __name__ == "__main__": - go() diff --git a/garf_models/test02_main4_compare.py b/garf_models/test02_main4_compare.py deleted file mode 100755 index b564355..0000000 --- a/garf_models/test02_main4_compare.py +++ /dev/null @@ -1,121 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -from pathlib import Path -from opengate.tests import utility -from garf.helpers import plot_spect_projection - - -def go(): - - folder_base = Path("test02") - folder_ref = folder_base / "reference" - folder_arf = folder_base / "arf" - folder_ff = folder_base / "free_flight" - stats_ref = utility.read_stat_file(folder_ref / "stats.txt") - - is_ok = True - scaling = 10 - - # ref vs ARF proj1 - print(f"Ref vs ARF (projection1) =======================") - im_name = "projection_1.mhd" - f1 = folder_ref / im_name - f2 = folder_arf / im_name - is_ok = ( - utility.assert_images( - f1, - f2, - stats_ref, - test_sad=False, - sum_tolerance=15, - scaleImageValuesFactor=scaling, - axis="x", - sad_profile_tolerance=20, - ) - and is_ok - ) - - o = folder_arf / "compare_1.pdf" - plt = plot_spect_projection(f1, f2, scaling=10, islice=117, wslice=3) - plt.savefig(o) - print(f"Profile figure : ", o) - - # ref vs ARF proj2 - print() - print(f"Ref vs ARF (projection2) =======================") - im_name = "projection_2.mhd" - f1 = folder_ref / im_name - f2 = folder_arf / im_name - is_ok = ( - utility.assert_images( - f1, - f2, - stats_ref, - test_sad=False, - sum_tolerance=15, - scaleImageValuesFactor=scaling, - axis="x", - sad_profile_tolerance=20, - ) - and is_ok - ) - - o = folder_arf / "compare_2.pdf" - plt = plot_spect_projection(f1, f2, scaling=10, islice=117, wslice=3) - plt.savefig(o) - print(f"Profile figure : ", o) - - # ref vs FF proj1 - print() - print(f"Ref vs FF (projection1) =======================") - im_name = "projection_1.mhd" - f1 = folder_ref / im_name - f2 = folder_ff / im_name - is_ok = ( - utility.assert_images( - f1, - f2, - stats_ref, - test_sad=False, - sum_tolerance=70, - scaleImageValuesFactor=scaling, - axis="x", - ) - and is_ok - ) - - o = folder_ff / "compare_1.pdf" - plt = plot_spect_projection(f1, f2, scaling=10, islice=117, wslice=3) - plt.savefig(o) - print(f"Profile figure : ", o) - - print() - print(f"Ref vs FF (projection1) =======================") - im_name = "projection_2.mhd" - f1 = folder_ref / im_name - f2 = folder_ff / im_name - is_ok = ( - utility.assert_images( - f1, - f2, - stats_ref, - test_sad=False, - sum_tolerance=70, - scaleImageValuesFactor=scaling, - axis="x", - ) - and is_ok - ) - - o = folder_ff / "compare_2.pdf" - plt = plot_spect_projection(f1, f2, scaling=10, islice=117, wslice=3) - plt.savefig(o) - print(f"Profile figure : ", o) - - utility.test_ok(is_ok) - - -# -------------------------------------------------------------------------- -if __name__ == "__main__": - go() From dea3f281e31c53818379ca80dcab4ee06f0ffeff Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Thu, 10 Jul 2025 08:51:45 +0200 Subject: [PATCH 08/20] Add MultiTask ResNet architecture with custom loss function Introduced `MultiTask_v3` model for dual-task predictions (acceptance and energy window classification) using a shared ResNet backbone. Added custom loss calculation in `MultiTask_v3_loss` to handle binary and categorical outputs with class balancing support. Refactored existing code for consistency and modularity. --- garf/garf_model.py | 147 ++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 145 insertions(+), 2 deletions(-) diff --git a/garf/garf_model.py b/garf/garf_model.py index 73900ce..90bf89a 100644 --- a/garf/garf_model.py +++ b/garf/garf_model.py @@ -2,6 +2,8 @@ import torch import torch.nn as nn +import numpy as np +import torch.nn.functional as F class Net_v1(nn.Module): @@ -23,7 +25,7 @@ class Net_v1(nn.Module): def __init__(self, H, L, n_ene_win): super(Net_v1, self).__init__() # Linear include Bias=True by default - self.fc1 = nn.Linear(3, H) + self.input_layer = nn.Linear(3, H) self.L = L self.fcts = nn.ModuleList() for i in range(L): @@ -31,10 +33,151 @@ def __init__(self, H, L, n_ene_win): self.fc3 = nn.Linear(H, n_ene_win) def forward(self, X): - X = self.fc1(X) # first layer + X = self.input_layer(X) # first layer X = torch.clamp(X, min=0) # relu for i in range(self.L): X = self.fcts[i](X) # hidden layers X = torch.clamp(X, min=0) # relu X = self.fc3(X) # output layer return X + + +class ResidualBlock(nn.Module): + """A simple residual block with two linear layers.""" + + def __init__(self, H): + super().__init__() + self.linear1 = nn.Linear(H, H) + self.linear2 = nn.Linear(H, H) + + def forward(self, x): + # Calculate the residual + residual = self.linear1(x) + residual = torch.clamp(residual, min=0) # ReLU + residual = self.linear2(residual) + + # Add the input to the residual (the "skip connection") + # and apply the final activation for this block + out = x + residual + out = torch.clamp(out, min=0) # ReLU + return out + + +class ResNet_v2(nn.Module): + """A ResNet-style architecture for the ARF problem.""" + + def __init__(self, H, L, n_ene_win): + super().__init__() + # Initial layer to project input from 3 dimensions to H dimensions + self.input_layer = nn.Linear(3, H) + + # A series of residual blocks + self.residual_layers = nn.ModuleList([ResidualBlock(H) for _ in range(L)]) + + # Final output layer + self.output_layer = nn.Linear(H, n_ene_win) + + def forward(self, x): + # Pass through the input layer and apply first activation + x = self.input_layer(x) + x = torch.clamp(x, min=0) # ReLU + + # Pass through all the residual blocks + for layer in self.residual_layers: + x = layer(x) + + # Final prediction + x = self.output_layer(x) + return x + + +class MultiTask_v3(nn.Module): + """ + A multi-task ResNet architecture for the ARF problem. + It has two output heads: one for acceptance and one for energy window classification. + """ + + def __init__(self, H, L, n_energy_windows): + super().__init__() + # Shared backbone + self.input_layer = nn.Linear(3, H) + self.residual_layers = nn.ModuleList([ResidualBlock(H) for _ in range(L)]) + + # Head 1: Predicts detection vs. non-detection (1 output logit) + self.acceptance_head = nn.Sequential( + nn.Linear(H, H // 2), nn.ReLU(), nn.Linear(H // 2, 1) + ) + + # Head 2: Predicts which detected window (n_energy_windows-1 outputs) + # We subtract 1 because we don't need to predict the "non-detected" class here. + self.energy_head = nn.Sequential( + nn.Linear(H, H // 2), nn.ReLU(), nn.Linear(H // 2, n_energy_windows - 1) + ) + + def forward(self, x): + # Pass through the shared backbone + x = self.input_layer(x) + x = torch.clamp(x, min=0) # ReLU + for layer in self.residual_layers: + x = layer(x) + + # Get predictions from each head + acceptance_logit = self.acceptance_head(x) + energy_logits = self.energy_head(x) + + return acceptance_logit, energy_logits + + +class MultiTask_v3_loss: + + def __init__(self, y_train, rr_factor, current_gpu_device): + print("init") + self.current_gpu_device = current_gpu_device + + # Get class weights for the acceptance loss (same hybrid method as before) + class_counts = np.bincount(y_train) + # Create binary weights: weight for non-detected (0) vs. detected (1) + weight_0 = 1.0 / (class_counts[0] + 1e-9) + weight_1 = 1.0 / (np.sum(class_counts[1:]) + 1e-9) + self.acceptance_weights = torch.tensor( + [weight_0, weight_1], dtype=torch.float + ).to(current_gpu_device) + if rr_factor > 1: + self.acceptance_weights[0] *= rr_factor + self.acceptance_weights /= torch.mean(self.acceptance_weights) + print(f"Acceptance loss weights: {self.acceptance_weights.cpu().numpy()}") + + def loss(self, Y_out, Y_true): + # Forward pass - model now returns two outputs + acceptance_logit, energy_logits = Y_out + + # --- Custom Loss Calculation --- + + # 1. Acceptance Loss (Binary) + # Create binary target: 0 if non-detected, 1 if detected + Y_binary_true = (Y_true > 0).float().view(-1, 1) + # Use BCEWithLogitsLoss which is numerically stable and takes class weights + # We need to compute pos_weight for the binary case + pos_weight = torch.tensor( + [self.acceptance_weights[1] / self.acceptance_weights[0]], + device=self.current_gpu_device, + ) + acceptance_loss_fn = torch.nn.BCEWithLogitsLoss(pos_weight=pos_weight) + loss1 = acceptance_loss_fn(acceptance_logit, Y_binary_true) + + # 2. Energy Loss (Categorical, only for detected events) + detected_mask = Y_true > 0 + if torch.sum(detected_mask) > 0: + # We subtract 1 from the labels because energy_head has (n-1) outputs + # e.g., window 1 -> class 0, window 2 -> class 1 + Y_energy_true = Y_true[detected_mask] - 1 + energy_logits_detected = energy_logits[detected_mask] + loss2 = F.cross_entropy(energy_logits_detected, Y_energy_true) + else: + # If no detected events in this batch, loss is zero + loss2 = 0.0 + + # 3. Total Loss + loss = loss1 + loss2 # You can weight these, e.g., loss1 + 0.5 * loss2 + + return loss From f3071a1b85ac3ddfccd2c94e2004f1394293a7a0 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Thu, 10 Jul 2025 08:51:55 +0200 Subject: [PATCH 09/20] Refactor and extend model compatibility in scatter prediction. Introduced support for multiple model types (e.g., MultiTask_v3, ResNet_v2, and xgboost) and enhanced predictions with modular functions. Added handling for plane axis initialization and improved geometric projection for detector plane interactions. Refactored functions for better abstraction and future extensibility. --- garf/garf_detector.py | 267 +++++++++++++++++++++++++++++++++++++++--- 1 file changed, 250 insertions(+), 17 deletions(-) diff --git a/garf/garf_detector.py b/garf/garf_detector.py index 006c35d..198bcea 100644 --- a/garf/garf_detector.py +++ b/garf/garf_detector.py @@ -1,9 +1,7 @@ # -*- coding: utf-8 -*- -import numpy as np -import torch from torch import Tensor import itk -from .garf_model import Net_v1 +from .garf_model import * from .helpers import get_gpu_device @@ -40,12 +38,23 @@ def load_nn(filename, verbose=True, gpu_mode="auto"): "Best epoch = {}".format(nn["optim"]["data"][best_epoch_eval]["epoch"]) ) + model_type = "Net_v1" + if "model_type" in model_data: + model_type = model_data["model_type"] + # prepare the model state = nn["optim"]["model_state"][best_epoch_eval] H = model_data["H"] n_ene_win = model_data["n_ene_win"] L = model_data["L"] - model = Net_v1(H, L, n_ene_win) + + if model_type == "Net_v1": + model = Net_v1(H, L, n_ene_win) + if model_type == "ResNet_v2": + model = ResNet_v2(H, L, n_ene_win) + if model_type == "MultiTask_v3": + model = MultiTask_v3(H, L, n_ene_win) + model.load_state_dict(state) return nn, model @@ -62,6 +71,7 @@ def __init__(self): self.initial_plane_rotation = None self.radius = None self.hit_slice_flag = False + self.plane_axis = None # computed self.plane_rotations = None @@ -139,6 +149,10 @@ def initialize(self, gantry_rotations): ) exit(-1) + # plane axis + if self.plane_axis is None: + self.plane_axis = [0, 1, 2] + # planes self.initialize_detector_plane_rotations(gantry_rotations) @@ -425,6 +439,7 @@ def image_from_coordinates_add_torch_hit_slice(self, img, vu): class GarfDetectorPlane: def __init__(self, garf_detector, center, rotation): self.garf_detector = garf_detector + self.plane_axis = self.garf_detector.plane_axis self.M = self.rotation_to_tensor(rotation) self.Mt = self.M.t() self.center = center @@ -486,7 +501,6 @@ def plane_intersection(self, batch): # two first coord pos_xy_rot = pos_xyz_rot[:, 0:2] - dir_xy_rot = dir_xyz_rot[:, 0:2] s = self.garf_detector.image_size_mm indexes_to_keep = torch.where( @@ -495,11 +509,13 @@ def plane_intersection(self, batch): )[0] # convert direction into theta/phi - # theta is acos(dy) - # phi is acos(dx) - nb = len(dir_xy_rot) - theta = torch.rad2deg(torch.arccos(dir_xy_rot[:, 1])).reshape((nb, 1)) - phi = torch.rad2deg(torch.arccos(dir_xy_rot[:, 0])).reshape((nb, 1)) + nb = len(dir_xyz_rot) + d_x_plane = dir_xyz_rot[:, self.plane_axis[0]] + d_y_plane = dir_xyz_rot[:, self.plane_axis[1]] + d_z_plane = dir_xyz_rot[:, self.plane_axis[2]] + theta = torch.rad2deg(torch.arccos(d_z_plane)).reshape((nb, 1)) + phi = torch.rad2deg(torch.arctan2(d_y_plane, d_x_plane)).reshape((nb, 1)) + angles = torch.concat((theta, phi), dim=1) batch = torch.concat( @@ -595,9 +611,15 @@ def nn_predict_numpy(model, model_data, x): # predict values vy_pred = model(vx) + # Apply correction to the logits BEFORE softmax + # Adding log(rr) to a logit is equivalent to multiplying its + # probability by rr after exponentiation. + if rr > 1: + vy_pred[:, 0] += np.log(rr) + # convert to numpy and normalize probabilities y_pred = normalize_logproba(vy_pred.data) - y_pred = normalize_proba_with_russian_roulette(y_pred, 0, rr) + # y_pred = normalize_proba_with_russian_roulette(y_pred, 0, rr) y_pred = y_pred.cpu().numpy() y_pred = y_pred.astype(np.float64) @@ -605,6 +627,74 @@ def nn_predict_numpy(model, model_data, x): return y_pred +def nn_predict_numpy_multitask(model, model_data, x): + """ + Apply the Multi-Task NN to predict y from x. + This version correctly handles the two-headed output. + """ + # Apply input normalization (same as before) + x_mean = model_data["x_mean"] + x_std = model_data["x_std"] + x = (x - x_mean) / x_std + + # Set device and model (same as before) + device = model_data.get("current_gpu_device", torch.device("cpu")) + model.to(device) + model.eval() # Set model to evaluation mode + + # Torch encapsulation + x = x.astype("float32") + vx = torch.from_numpy(x).to(device) + + # --- New Multi-Task Prediction Logic --- + + with torch.no_grad(): # Disable gradient calculation for inference + # 1. Get the two outputs from the model + acceptance_logit, energy_logits = model(vx) + + # 2. Calculate probability of detection using the sigmoid function + # This is P(detection) + p_acceptance = torch.sigmoid(acceptance_logit) + + # 3. Calculate conditional probabilities for each energy window using softmax + # This is P(window k | detected) + p_energy_windows = torch.softmax(energy_logits, dim=1) + + # 4. Combine the probabilities to get the final prediction + p_acceptance = p_acceptance.cpu().numpy() + p_energy_windows = p_energy_windows.cpu().numpy() + + # The number of final windows is 1 (for non-detected) + number of energy heads + n_samples = x.shape[0] + n_total_windows = p_energy_windows.shape[1] + 1 + y_pred = np.zeros((n_samples, n_total_windows), dtype=np.float64) + + # Probability of non-detection (window 0) is 1 - P(detection) + y_pred[:, 0] = 1.0 - p_acceptance.flatten() + + # Probability of a given detected window k is P(detection) * P(window k | detected) + for k in range(p_energy_windows.shape[1]): + y_pred[:, k + 1] = p_acceptance.flatten() * p_energy_windows[:, k] + + return y_pred + + +def xgb_predict_numpy(model, model_data, x): + """ + Apply a trained XGBoost model to predict y from x. + """ + # Apply the same normalization used during training + x_mean = model_data["x_mean"] + x_std = model_data["x_std"] + x_normalized = (x - x_mean) / x_std + + # Predict probabilities directly; no RR correction needed + # as the model was trained with weighted loss. + y_pred = model.predict_proba(x_normalized) + + return y_pred.astype(np.float64) + + def compute_angle_offset_torch(angles, length): """ compute the x,y offset according to the angle @@ -711,7 +801,7 @@ def image_from_coordinates_add_numpy(img, u, v, w_pred, hit_slice=False): img[i, uv16Bins[chx > tiny, 0], uv16Bins[chx > tiny, 1]] += chx[chx > tiny] -def arf_plane_intersection(batch, plane, image_plane_size_mm): +def arf_plane_intersection(batch, plane, image_plane_size_mm, plane_axis): """ Project the x points (Ekine X Y Z dX dY dZ) on the image plane defined by plane_U, plane_V, plane_center, plane_normal @@ -792,14 +882,18 @@ def arf_plane_intersection(batch, plane, image_plane_size_mm): # two first coord of dir dx = dir_xyz_rot[:, 0] dy = dir_xyz_rot[:, 1] + dz = dir_xyz_rot[:, 2] # FIXME -> clip arcos -1;1 ? # convert direction into theta/phi - # theta is acos(dy) - # phi is acos(dx) - theta = np.degrees(np.arccos(dy)).reshape((nb, 1)) - phi = np.degrees(np.arccos(dx)).reshape((nb, 1)) + dirs = np.stack((dx, dy, dz), axis=-1) + d_x_plane = dirs[:, plane_axis[0]] + d_y_plane = dirs[:, plane_axis[1]] + d_z_plane = dirs[:, plane_axis[2]] + theta = np.degrees(np.arccos(d_z_plane)).reshape((nb, 1)) + phi = np.degrees(np.arctan2(d_y_plane, d_x_plane)).reshape((nb, 1)) + y = np.concatenate((y, theta), axis=1) y = np.concatenate((y, phi), axis=1) @@ -863,7 +957,7 @@ def arf_from_points_to_image_counts_OLD( return u, v, w_pred -def arf_from_points_to_image_counts( +def arf_from_points_to_image_counts_OLD2( projected_batch, # 5D: 2 plane coordinates, 2 angles, 1 energy, 1 weight model, # ARF neural network model model_data, # associated model data @@ -920,3 +1014,142 @@ def arf_from_points_to_image_counts( ) return u, v, w_pred + + +def arf_from_points_to_image_counts_OLD3( + projected_batch, # Now 6D/7D: pos_x, pos_y, dir_x, dir_y, dir_z, E, [weight] + model, + model_data, + distance_to_crystal, + image_plane_size_mm, + image_plane_size_pixel, + image_plane_spacing, +): + """ + Input: position, direction on the detector plane, energy. + This version performs a correct geometric projection. + """ + + # --- 1. Predict scatter probabilities with the NN --- + + # Get directions (dx, dy, dz) and energy for NN input + dirs = projected_batch[:, 2:5] + energy = projected_batch[:, 5:6] + + # Calculate angles internally for the NN + theta = np.degrees(np.arccos(dirs[:, 2])) + phi = np.degrees(np.arctan2(dirs[:, 1], dirs[:, 0])) + + # Assemble NN input (theta, phi, E) and predict + ax = np.column_stack((theta, phi, energy)) + w_pred = nn_predict_numpy(model, model_data, ax) + + # Apply particle weights if they exist + if projected_batch.shape[1] == 7: + weights = projected_batch[:, 6] + w_pred = w_pred * weights[:, np.newaxis] + + # --- 2. Calculate final position on the crystal --- + + # Get initial position (px, py) + cx = projected_batch[:, 0:2] + + # Correctly project the trajectory over distance_to_crystal + dir_z = dirs[:, 2] + # Handle particles traveling parallel to the plane to avoid division by zero + mask = np.abs(dir_z) > 1e-9 + offset = np.zeros_like(cx) + offset[mask, 0] = distance_to_crystal * (dirs[mask, 0] / dir_z[mask]) # d * dx/dz + offset[mask, 1] = distance_to_crystal * (dirs[mask, 1] / dir_z[mask]) # d * dy/dz + + # The final position on the crystal plane + final_pos = cx + offset + + # --- 3. Convert coordinates to image pixels --- + + # Convert mm coordinates to pixel indices + coord = ( + final_pos + image_plane_size_mm / 2 - image_plane_spacing / 2 + ) / image_plane_spacing + coord = np.around(coord).astype(int) + + # Separate into u, v coordinates + v = coord[:, 0] + u = coord[:, 1] + + # Remove points that fall outside the image dimensions + u, v, w_pred = remove_out_of_image_boundaries_numpy( + u, v, w_pred, image_plane_size_pixel + ) + + return u, v, w_pred + + +def arf_from_points_to_image_counts( + projected_batch, + model, + model_data, + distance_to_crystal, + image_plane_size_mm, + image_plane_size_pixel, + image_plane_spacing, +): + """ + Input: position, direction on the detector plane, energy. + This version is model-agnostic (PyTorch or XGBoost). + """ + + # --- 1. Predict scatter probabilities with the correct model --- + + # Get directions and energy for model input + dirs = projected_batch[:, 2:5] + energy = projected_batch[:, 5:6] + + # Calculate angles from directions + theta = np.degrees(np.arccos(np.clip(dirs[:, 2], -1, 1))) + phi = np.degrees(np.arctan2(dirs[:, 1], dirs[:, 0])) + + # Assemble input (theta, phi, E) + ax = np.column_stack((theta, phi, energy)) + + # Conditionally call the correct prediction function + model_type = model_data.get("model_type", "Net_v1") + if model_type == "xgboost": + w_pred = xgb_predict_numpy(model, model_data, ax) + elif model_type == "MultiTask_v3": + w_pred = nn_predict_numpy_multitask(model, model_data, ax) + else: + w_pred = nn_predict_numpy(model, model_data, ax) + + # Apply particle weights from the simulation if they exist + if projected_batch.shape[1] == 7: + weights = projected_batch[:, 6] + w_pred = w_pred * weights[:, np.newaxis] + + # --- The rest of the function (geometric projection) remains the same --- + + # Get initial position (px, py) + cx = projected_batch[:, 0:2] + + # Correctly project the trajectory over distance_to_crystal + dir_z = dirs[:, 2] + mask = np.abs(dir_z) > 1e-9 + offset = np.zeros_like(cx) + offset[mask, 0] = distance_to_crystal * (dirs[mask, 0] / dir_z[mask]) + offset[mask, 1] = distance_to_crystal * (dirs[mask, 1] / dir_z[mask]) + final_pos = cx + offset + + # Convert mm coordinates to pixel indices + coord = ( + final_pos + image_plane_size_mm / 2 - image_plane_spacing / 2 + ) / image_plane_spacing + coord = np.around(coord).astype(int) + + # Separate into u, v coordinates and remove out-of-bounds points + v = coord[:, 0] + u = coord[:, 1] + u, v, w_pred = remove_out_of_image_boundaries_numpy( + u, v, w_pred, image_plane_size_pixel + ) + + return u, v, w_pred From 182c5390eb90fa036b2cff22720ad35cd91dc466 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Thu, 10 Jul 2025 08:52:17 +0200 Subject: [PATCH 10/20] Remove unused imports from garf_plot_training_dataset.py Eliminated unnecessary modules like sys, cm, np, uproot, and ntpath to streamline the script. This improves code readability and reduces potential overhead from unused dependencies. --- garf/bin/garf_plot_training_dataset.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/garf/bin/garf_plot_training_dataset.py b/garf/bin/garf_plot_training_dataset.py index 0f0914d..ea6c109 100755 --- a/garf/bin/garf_plot_training_dataset.py +++ b/garf/bin/garf_plot_training_dataset.py @@ -1,13 +1,8 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- -import sys import garf import matplotlib.pyplot as plt -from matplotlib import cm -import numpy as np -import uproot -import ntpath import click # ----------------------------------------------------------------------------- From 8fda77767679d6ee75ab581c69271a32fdea408d Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Thu, 10 Jul 2025 08:52:29 +0200 Subject: [PATCH 11/20] Add rr and epoch options to garf_train CLI This update introduces two new CLI options, `--rr` and `--epoch`, allowing users to override the corresponding parameters in the config file. These additions enhance training flexibility and enable quick parameter adjustments directly on the command line. --- garf/bin/garf_train.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/garf/bin/garf_train.py b/garf/bin/garf_train.py index 86648c8..00b85d4 100755 --- a/garf/bin/garf_train.py +++ b/garf/bin/garf_train.py @@ -16,7 +16,9 @@ @click.argument("data") @click.argument("output") @click.option("--progress-bar/--no-progress-bar", default=True) -def garf_train(param, data, output, progress_bar): +@click.option("--rr", default=None, help="RR value (overwrite the one in the param file)") +@click.option("--epoch", default=None, help="Nb of epoch (overwrite the one in the param file)") +def garf_train(param, data, output, rr, epoch, progress_bar): """ \b Train a ARF-nn (neural network) from a training dataset. @@ -34,6 +36,14 @@ def garf_train(param, data, output, progress_bar): params = json.loads(param_file) params["progress_bar"] = progress_bar + # RR ? + if rr is not None: + params["RR"] = float(rr) + + # epoch ? + if epoch is not None: + params["epoch_max"] = int(epoch) + # Print info print("Training dataset", data_filename) garf.print_training_dataset_info(data, params["RR"]) From 5aab539008d31557316d52582cf02f581b9828d3 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Thu, 10 Jul 2025 08:53:53 +0200 Subject: [PATCH 12/20] Remove default value for 'rr' parameter in print_training_dataset_info The function no longer sets a default value for the 'rr' parameter, requiring it to be explicitly provided. This change ensures greater clarity and avoids unintended behavior from hidden defaults. --- garf/helpers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/garf/helpers.py b/garf/helpers.py index 8a48d44..6b45252 100644 --- a/garf/helpers.py +++ b/garf/helpers.py @@ -45,7 +45,7 @@ def load_training_dataset(filename): return data, theta, phi, E, w -def print_training_dataset_info(data, rr=40): +def print_training_dataset_info(data, rr): """ Print training dataset information """ From 40ed918ac38b2b5b4914680d4cd260872a6de93b Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Thu, 10 Jul 2025 08:54:00 +0200 Subject: [PATCH 13/20] Refactor training logic and enhance model flexibility Replaced hardcoded model configuration with dynamic type selection, supporting multiple architectures (Net_v1, ResNet_v2, MultiTask_v3). Optimized DataLoader setup, adjusted training loop parameters, and introduced early stopping in a new `train_nn_TEST` function for efficient training workflows. --- garf/garf_train.py | 143 ++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 130 insertions(+), 13 deletions(-) diff --git a/garf/garf_train.py b/garf/garf_train.py index 50b6f5a..85999ed 100644 --- a/garf/garf_train.py +++ b/garf/garf_train.py @@ -6,7 +6,7 @@ from torch import Tensor import copy from tqdm import tqdm -from .garf_model import Net_v1 +from .garf_model import * from .helpers import get_gpu_device @@ -297,13 +297,11 @@ def train_nn(x_train, y_train, params): print(f"Device GPU type is {current_gpu_mode}") # Batch parameters - batch_per_epoch = model_data["batch_per_epoch"] batch_size = model_data["batch_size"] epoch_store_every = model_data["epoch_store_every"] # DataLoader print("Data loader batch_size", batch_size) - print("Data loader batch_per_epoch", batch_per_epoch) train_data2 = np.column_stack((x_train, y_train)) if current_gpu_mode == "mps": print("With device mps (gpu), convert data to float32", train_data2.dtype) @@ -312,8 +310,8 @@ def train_nn(x_train, y_train, params): train_loader2 = DataLoader( train_data2, batch_size=batch_size, - num_workers=2, - # pin_memory=True, + num_workers=8, + pin_memory=True, # shuffle=True, # if false ~20% faster, seems identical shuffle=False, # if false ~20% faster, seems identical drop_last=True, @@ -322,7 +320,20 @@ def train_nn(x_train, y_train, params): # Create the main NN H = model_data["H"] L = model_data["L"] - model = Net_v1(H, L, n_ene_win) + rr_factor = params["RR"] + + model_type = "Net_v1" + loss_function = F.cross_entropy + if "model_type" in model_data: + model_type = model_data["model_type"] + if model_type == "Net_v1": + model = Net_v1(H, L, n_ene_win) + if model_type == "ResNet_v2": + model = ResNet_v2(H, L, n_ene_win) + if model_type == "MultiTask_v3": + model = MultiTask_v3(H, L, n_ene_win) + l = MultiTask_v3_loss(y_train, rr_factor, current_gpu_device) + loss_function = l.loss # Create the optimizer # optimizer, scheduler = nn_get_optimiser(model_data, model) @@ -331,7 +342,7 @@ def train_nn(x_train, y_train, params): optimizer = torch.optim.Adam( model.parameters(), lr=learning_rate, - # weight_decay=1e-4, + weight_decay=1e-4, ) # decreasing learning_rate scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau( @@ -376,7 +387,7 @@ def train_nn(x_train, y_train, params): for batch_idx, data in enumerate(train_loader2): x = data[:, 0:3] y = data[:, 3] - X = Tensor(x.to(model.fc1.weight.dtype)).to(current_gpu_device) + X = Tensor(x.to(model.input_layer.weight.dtype)).to(current_gpu_device) Y = Tensor(y).to(current_gpu_device).long() # Forward pass @@ -384,7 +395,7 @@ def train_nn(x_train, y_train, params): # Compute expected loss # combines log_softmax and nll_loss in a single function - loss = F.cross_entropy(Y_out, Y) + loss = loss_function(Y_out, Y) # Backward pass loss.backward() @@ -396,9 +407,6 @@ def train_nn(x_train, y_train, params): train_loss += loss.data.item() * batch_size n_samples_processed += batch_size - # Stop when batch_per_epoch is reach - if batch_idx == params["batch_per_epoch"]: - break # end for loop train_loader # end of train @@ -435,7 +443,7 @@ def train_nn(x_train, y_train, params): print("Store weights", epoch) optim_data["epoch"] = epoch optim_data["train_loss"] = train_loss - state = copy.deepcopy(model.state_dict()) + state = model.state_dict() nn["optim"]["model_state"].append(state) nn["optim"]["data"].append(optim_data) best_epoch_index = epoch @@ -455,3 +463,112 @@ def train_nn(x_train, y_train, params): model_data["best_loss"] = best_loss return nn + + +def train_nn_TEST(x_train, y_train, params): + """ + An optimized function to train the ARF neural network. + """ + # 1. PREPARE DATA AND MODEL PARAMETERS + # =================================================================== + print("Preparing data and model parameters...") + x_train, y_train, model_data, N = nn_prepare_data(x_train, y_train, params) + y_vals, y_train = np.unique(y_train, return_inverse=True) + model_data["n_ene_win"] = len(y_vals) + current_gpu_mode, current_gpu_device = get_gpu_device(params["gpu_mode"]) + model_data["current_gpu_mode"] = current_gpu_mode + print_nn_params(model_data) + + # 2. CREATE EFFICIENT DATALOADER + # =================================================================== + # This is the most efficient method for small, in-memory datasets. + print("Creating PyTorch TensorDataset and DataLoader...") + + # Convert NumPy arrays to PyTorch Tensors ONCE. + x_tensor = torch.from_numpy(x_train.astype(np.float32)) + y_tensor = torch.from_numpy(y_train.astype(np.int64)) + + # Create a TensorDataset + train_dataset = torch.utils.data.TensorDataset(x_tensor, y_tensor) + + # Use DataLoader with num_workers=0 (no multiprocessing overhead) + # and pin_memory=True (for faster CPU to CUDA transfer). + train_loader = DataLoader( + dataset=train_dataset, + batch_size=model_data["batch_size"], + num_workers=0, + pin_memory=True, + shuffle=True, # Shuffle is good practice for training + drop_last=True, + ) + + # 3. SETUP MODEL, OPTIMIZER, and LOSS + # =================================================================== + print("Setting up model, optimizer, and loss function...") + model = Net_v1(model_data["H"], model_data["L"], model_data["n_ene_win"]) + model.to(current_gpu_device) + + optimizer = torch.optim.Adam(model.parameters(), lr=model_data["learning_rate"]) + scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(optimizer, "min", patience=5) + + # (Optional, but recommended) - Use the weighted loss for better results + # class_counts = np.bincount(y_train) + # class_weights = 1.0 / (class_counts + 1e-9) + # class_weights[0] *= model_data["RR"] + # class_weights = class_weights / np.mean(class_weights) + # class_weights = torch.tensor(class_weights, dtype=torch.float).to(current_gpu_device) + + # 4. MAIN TRAINING LOOP + # =================================================================== + print("\nStarting optimized training...") + nn = {"model_data": model_data, "optim": {"model_state": [], "data": []}} + best_loss = np.Inf + best_epoch = -1 + epochs_no_improve = 0 + pbar = tqdm(range(model_data["epoch_max"]), disable=not params["progress_bar"]) + for epoch in pbar: + model.train() + total_loss = 0 + + # An epoch is now one full pass over the entire dataset + for x_batch, y_batch in train_loader: + # Move data to GPU + X = x_batch.to(current_gpu_device) + Y = y_batch.to(current_gpu_device) + + # Standard forward/backward pass + optimizer.zero_grad() + Y_out = model(X) + + # Use unweighted loss as requested (or uncomment weighted version) + loss = F.cross_entropy(Y_out, Y) + # loss = F.cross_entropy(Y_out, Y, weight=class_weights) + + loss.backward() + optimizer.step() + total_loss += loss.item() + + # After each epoch, calculate average loss + avg_loss = total_loss / len(train_loader) + scheduler.step(avg_loss) + pbar.set_description(f"Epoch {epoch + 1}, Loss: {avg_loss:.5f}") + + # Check for improvement (for early stopping) + if avg_loss < best_loss: + best_loss = avg_loss + best_epoch = epoch + epochs_no_improve = 0 + # Store the best model state efficiently + nn["optim"]["model_state"] = [model.state_dict()] + nn["optim"]["data"] = [{"epoch": epoch, "train_loss": avg_loss}] + nn["model_data"]["best_epoch"] = epoch + nn["model_data"]["best_loss"] = best_loss + else: + epochs_no_improve += 1 + + if epochs_no_improve >= model_data["early_stopping"]: + print(f"\nEarly stopping triggered after {epoch + 1} epochs.") + break + + print(f"\nTraining done. Best loss = {best_loss:.5f} at epoch {best_epoch + 1}") + return nn From 844b56bf920f6fda55f94df10b2d80f39a479dca Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Mon, 28 Jul 2025 16:59:55 +0200 Subject: [PATCH 14/20] Refactor: Rename `input_layer` to `fc1` for consistency Renamed the `input_layer` variable to `fc1` across all relevant classes and methods for improved naming consistency with other layers. This change enhances code readability and aligns variable naming conventions in the model definitions. --- garf/garf_model.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/garf/garf_model.py b/garf/garf_model.py index 90bf89a..ea9c95b 100644 --- a/garf/garf_model.py +++ b/garf/garf_model.py @@ -25,7 +25,7 @@ class Net_v1(nn.Module): def __init__(self, H, L, n_ene_win): super(Net_v1, self).__init__() # Linear include Bias=True by default - self.input_layer = nn.Linear(3, H) + self.fc1 = nn.Linear(3, H) self.L = L self.fcts = nn.ModuleList() for i in range(L): @@ -33,7 +33,7 @@ def __init__(self, H, L, n_ene_win): self.fc3 = nn.Linear(H, n_ene_win) def forward(self, X): - X = self.input_layer(X) # first layer + X = self.fc1(X) # first layer X = torch.clamp(X, min=0) # relu for i in range(self.L): X = self.fcts[i](X) # hidden layers @@ -69,7 +69,7 @@ class ResNet_v2(nn.Module): def __init__(self, H, L, n_ene_win): super().__init__() # Initial layer to project input from 3 dimensions to H dimensions - self.input_layer = nn.Linear(3, H) + self.fc1 = nn.Linear(3, H) # A series of residual blocks self.residual_layers = nn.ModuleList([ResidualBlock(H) for _ in range(L)]) @@ -79,7 +79,7 @@ def __init__(self, H, L, n_ene_win): def forward(self, x): # Pass through the input layer and apply first activation - x = self.input_layer(x) + x = self.fc1(x) x = torch.clamp(x, min=0) # ReLU # Pass through all the residual blocks @@ -100,7 +100,7 @@ class MultiTask_v3(nn.Module): def __init__(self, H, L, n_energy_windows): super().__init__() # Shared backbone - self.input_layer = nn.Linear(3, H) + self.fc1 = nn.Linear(3, H) self.residual_layers = nn.ModuleList([ResidualBlock(H) for _ in range(L)]) # Head 1: Predicts detection vs. non-detection (1 output logit) @@ -116,7 +116,7 @@ def __init__(self, H, L, n_energy_windows): def forward(self, x): # Pass through the shared backbone - x = self.input_layer(x) + x = self.fc1(x) x = torch.clamp(x, min=0) # ReLU for layer in self.residual_layers: x = layer(x) From 8b18dd998512005ba703f4006a9662bca3b413ca Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Mon, 28 Jul 2025 17:00:03 +0200 Subject: [PATCH 15/20] Update angle parametrization and fix tensor type reference Introduced a new angle parametrization using 'atan2' for better computation clarity. Also corrected the tensor type reference in the training loop to use `fc1.weight.dtype` instead of `input_layer.weight.dtype`. --- garf/garf_train.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/garf/garf_train.py b/garf/garf_train.py index 85999ed..f8336fe 100644 --- a/garf/garf_train.py +++ b/garf/garf_train.py @@ -49,6 +49,10 @@ def nn_prepare_data(x_train, y_train, params): model_data["x_std"] = x_std model_data["N"] = N + # this flag indicate that we use the new version of angle parametrisation + # (acos + atan2 instead of acos + acos) + model_data["angle_param"] = 'atan2' + # copy param except comments for i in params: if not i[0] == "#": @@ -387,7 +391,7 @@ def train_nn(x_train, y_train, params): for batch_idx, data in enumerate(train_loader2): x = data[:, 0:3] y = data[:, 3] - X = Tensor(x.to(model.input_layer.weight.dtype)).to(current_gpu_device) + X = Tensor(x.to(model.fc1.weight.dtype)).to(current_gpu_device) Y = Tensor(y).to(current_gpu_device).long() # Forward pass From 8b6637740502dae486b8237d043ab9b3d2e12728 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Mon, 28 Jul 2025 17:00:25 +0200 Subject: [PATCH 16/20] Refactor angle parametrization logic and remove old functions Introduced configurable angle parameterization (acos/atan2) and updated logic accordingly. Deprecated and removed legacy functions (`arf_from_points_to_image_counts_OLD*`) to streamline the codebase. Refactored function and method names to ensure consistency and clarity in angle computation workflows. --- garf/garf_detector.py | 225 +++++------------------------------------- 1 file changed, 22 insertions(+), 203 deletions(-) diff --git a/garf/garf_detector.py b/garf/garf_detector.py index 198bcea..607bd91 100644 --- a/garf/garf_detector.py +++ b/garf/garf_detector.py @@ -38,10 +38,18 @@ def load_nn(filename, verbose=True, gpu_mode="auto"): "Best epoch = {}".format(nn["optim"]["data"][best_epoch_eval]["epoch"]) ) + # which net type ? model_type = "Net_v1" if "model_type" in model_data: model_type = model_data["model_type"] + # which angle parametrisation ? + # old one = acos acos + # new one = acos atan2 + angle_param = 'acos' + if 'angle_param' not in model_data: + model_data['angle_param'] = angle_param + # prepare the model state = nn["optim"]["model_state"][best_epoch_eval] H = model_data["H"] @@ -264,7 +272,7 @@ def arf_plane_init_numpy(self, rotation, nb): def project_to_planes_torch(self, batch, projected_points): for i, detector_plane in enumerate(self.detector_planes): - projected_batch = detector_plane.plane_intersection(batch) + projected_batch = detector_plane.plane_intersection_torch(batch) self.build_image_from_projected_points_torch(projected_batch, i) def project_to_planes_numpy(self, batch, i, planes, projected_points, data_img): @@ -459,7 +467,7 @@ def rotation_to_tensor(self, m): t = t.to(self.garf_detector.current_gpu_device) return t - def plane_intersection(self, batch): + def plane_intersection_torch(self, batch): # See arf_plane_intersection # get energy, position and direction @@ -516,6 +524,10 @@ def plane_intersection(self, batch): theta = torch.rad2deg(torch.arccos(d_z_plane)).reshape((nb, 1)) phi = torch.rad2deg(torch.arctan2(d_y_plane, d_x_plane)).reshape((nb, 1)) + # FIXME previous angle parametrisation ? + #theta = torch.rad2deg(torch.arccos(dir_xy_rot[:, 1])).reshape((nb, 1)) + #phi = torch.rad2deg(torch.arccos(dir_xy_rot[:, 0])).reshape((nb, 1)) + angles = torch.concat((theta, phi), dim=1) batch = torch.concat( @@ -704,6 +716,8 @@ def compute_angle_offset_torch(angles, length): cos_theta = torch.cos(angles_rad[:, 0]) cos_phi = torch.cos(angles_rad[:, 1]) + print('FIXME acos !! ') + # see in Gate_NN_ARF_Actor, line "phi = acos(dir.x())/degree;" tx = length * cos_phi # see in Gate_NN_ARF_Actor, line "theta = acos(dir.y())/degree;" @@ -713,23 +727,6 @@ def compute_angle_offset_torch(angles, length): return t -def compute_angle_offset_numpy(angles, length): - """ - compute the x,y offset according to the angle - """ - angles_rad = np.deg2rad(angles) - cos_theta = np.cos(angles_rad[:, 0]) - cos_phi = np.cos(angles_rad[:, 1]) - - # see in Gate_NN_ARF_Actor, line "phi = acos(dir.x())/degree;" - tx = length * cos_phi - # see in Gate_NN_ARF_Actor, line "theta = acos(dir.y())/degree;" - ty = length * cos_theta - t = np.column_stack((tx, ty)) - - return t - - def normalize_proba_with_russian_roulette(w_pred, channel, rr): """ Consider rr times the values for the energy windows channel @@ -903,188 +900,6 @@ def arf_plane_intersection(batch, plane, image_plane_size_mm, plane_axis): return batch -def arf_from_points_to_image_counts_OLD( - projected_batch, # 5D: 2 plane coordinates, 2 angles, 1 energy - model, # ARF neural network model - model_data, # associated model data - distance_to_crystal, # from detection plane to crystal center - image_plane_size_mm, # image plane in mm - image_plane_size_pixel, # image plane in pixel - image_plane_spacing, -): # image plane spacing - """ - Input : position, direction on the detector plane, energy - Compute - - garf.nn_predict - - garf.compute_angle_offset - - garf.remove_out_of_image_boundaries2 - - Used in 1) GarfDetector class and 2) gate ARFActor - - """ - - # get the two angles and the energy - ax = projected_batch[:, 2:5] - - # predict weights - w_pred = nn_predict_numpy(model, model_data, ax) - - # Get the two first columns = points coordinates - cx = projected_batch[:, 0:2] - - # Get the two next columns = angles - angles = projected_batch[:, 2:4] - - # Take angle into account: consider position at collimator + half crystal - t = compute_angle_offset_numpy(angles, distance_to_crystal) - cx = cx + t - - # convert coord to pixel - coord = ( - cx + image_plane_size_mm / 2 - image_plane_spacing / 2 - ) / image_plane_spacing - coord = np.around(coord).astype(int) - - # why vu and not uv ? - v = coord[:, 0] - u = coord[:, 1] - - # remove points outside the image - u, v, w_pred = remove_out_of_image_boundaries_numpy( - u, v, w_pred, image_plane_size_pixel - ) - - return u, v, w_pred - - -def arf_from_points_to_image_counts_OLD2( - projected_batch, # 5D: 2 plane coordinates, 2 angles, 1 energy, 1 weight - model, # ARF neural network model - model_data, # associated model data - distance_to_crystal, # from detection plane to crystal center - image_plane_size_mm, # image plane in mm - image_plane_size_pixel, # image plane in pixel - image_plane_spacing, -): # image plane spacing - """ - Input : position, direction on the detector plane, energy - Compute - - garf.nn_predict - - garf.compute_angle_offset - - garf.remove_out_of_image_boundaries2 - - Used in 1) GarfDetector class and 2) gate ARFActor - - """ - - # get the two angles and the energy - ax = projected_batch[:, 2:5] - - # predict weights - w_pred = nn_predict_numpy(model, model_data, ax) - - # particle weight ? - if projected_batch.shape[1] == 6: - weights = projected_batch[:, 5] - w_pred = w_pred * weights[:, np.newaxis] - - # Get the two first columns = points coordinates - cx = projected_batch[:, 0:2] - - # Get the two next columns = angles - angles = projected_batch[:, 2:4] - - # Take angle into account: consider position at collimator + half crystal - t = compute_angle_offset_numpy(angles, distance_to_crystal) - cx = cx + t - - # convert coord to pixel - coord = ( - cx + image_plane_size_mm / 2 - image_plane_spacing / 2 - ) / image_plane_spacing - coord = np.around(coord).astype(int) - - # why vu and not uv ? - v = coord[:, 0] - u = coord[:, 1] - - # remove points outside the image - u, v, w_pred = remove_out_of_image_boundaries_numpy( - u, v, w_pred, image_plane_size_pixel - ) - - return u, v, w_pred - - -def arf_from_points_to_image_counts_OLD3( - projected_batch, # Now 6D/7D: pos_x, pos_y, dir_x, dir_y, dir_z, E, [weight] - model, - model_data, - distance_to_crystal, - image_plane_size_mm, - image_plane_size_pixel, - image_plane_spacing, -): - """ - Input: position, direction on the detector plane, energy. - This version performs a correct geometric projection. - """ - - # --- 1. Predict scatter probabilities with the NN --- - - # Get directions (dx, dy, dz) and energy for NN input - dirs = projected_batch[:, 2:5] - energy = projected_batch[:, 5:6] - - # Calculate angles internally for the NN - theta = np.degrees(np.arccos(dirs[:, 2])) - phi = np.degrees(np.arctan2(dirs[:, 1], dirs[:, 0])) - - # Assemble NN input (theta, phi, E) and predict - ax = np.column_stack((theta, phi, energy)) - w_pred = nn_predict_numpy(model, model_data, ax) - - # Apply particle weights if they exist - if projected_batch.shape[1] == 7: - weights = projected_batch[:, 6] - w_pred = w_pred * weights[:, np.newaxis] - - # --- 2. Calculate final position on the crystal --- - - # Get initial position (px, py) - cx = projected_batch[:, 0:2] - - # Correctly project the trajectory over distance_to_crystal - dir_z = dirs[:, 2] - # Handle particles traveling parallel to the plane to avoid division by zero - mask = np.abs(dir_z) > 1e-9 - offset = np.zeros_like(cx) - offset[mask, 0] = distance_to_crystal * (dirs[mask, 0] / dir_z[mask]) # d * dx/dz - offset[mask, 1] = distance_to_crystal * (dirs[mask, 1] / dir_z[mask]) # d * dy/dz - - # The final position on the crystal plane - final_pos = cx + offset - - # --- 3. Convert coordinates to image pixels --- - - # Convert mm coordinates to pixel indices - coord = ( - final_pos + image_plane_size_mm / 2 - image_plane_spacing / 2 - ) / image_plane_spacing - coord = np.around(coord).astype(int) - - # Separate into u, v coordinates - v = coord[:, 0] - u = coord[:, 1] - - # Remove points that fall outside the image dimensions - u, v, w_pred = remove_out_of_image_boundaries_numpy( - u, v, w_pred, image_plane_size_pixel - ) - - return u, v, w_pred - - def arf_from_points_to_image_counts( projected_batch, model, @@ -1106,8 +921,12 @@ def arf_from_points_to_image_counts( energy = projected_batch[:, 5:6] # Calculate angles from directions - theta = np.degrees(np.arccos(np.clip(dirs[:, 2], -1, 1))) - phi = np.degrees(np.arctan2(dirs[:, 1], dirs[:, 0])) + if model_data['angle_param'] == 'atan2': + theta = np.degrees(np.arccos(np.clip(dirs[:, 2], -1, 1))) + phi = np.degrees(np.arctan2(dirs[:, 1], dirs[:, 0])) + else: + theta = np.degrees(np.arccos(dirs[:, 1])) + phi = np.degrees(np.arccos(dirs[:, 0])) # Assemble input (theta, phi, E) ax = np.column_stack((theta, phi, energy)) From 3a028e5bce8bf03f456403f9603ec9c483b038d0 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Tue, 29 Jul 2025 08:07:25 +0200 Subject: [PATCH 17/20] Standardize string quotes and clean up comments. Replaced single quotes with double quotes for consistency. Removed outdated or redundant comments and unnecessary line breaks to improve code readability and maintainability. --- garf/garf_detector.py | 23 ++++++++--------------- 1 file changed, 8 insertions(+), 15 deletions(-) diff --git a/garf/garf_detector.py b/garf/garf_detector.py index 607bd91..0818da2 100644 --- a/garf/garf_detector.py +++ b/garf/garf_detector.py @@ -46,9 +46,9 @@ def load_nn(filename, verbose=True, gpu_mode="auto"): # which angle parametrisation ? # old one = acos acos # new one = acos atan2 - angle_param = 'acos' - if 'angle_param' not in model_data: - model_data['angle_param'] = angle_param + angle_param = "acos" + if "angle_param" not in model_data: + model_data["angle_param"] = angle_param # prepare the model state = nn["optim"]["model_state"][best_epoch_eval] @@ -525,8 +525,8 @@ def plane_intersection_torch(self, batch): phi = torch.rad2deg(torch.arctan2(d_y_plane, d_x_plane)).reshape((nb, 1)) # FIXME previous angle parametrisation ? - #theta = torch.rad2deg(torch.arccos(dir_xy_rot[:, 1])).reshape((nb, 1)) - #phi = torch.rad2deg(torch.arccos(dir_xy_rot[:, 0])).reshape((nb, 1)) + # theta = torch.rad2deg(torch.arccos(dir_xy_rot[:, 1])).reshape((nb, 1)) + # phi = torch.rad2deg(torch.arccos(dir_xy_rot[:, 0])).reshape((nb, 1)) angles = torch.concat((theta, phi), dim=1) @@ -542,9 +542,6 @@ def plane_intersection_torch(self, batch): return batch -# noinspection PyUnreachableCode - - def normalize_logproba(x): """ Convert un-normalized log probabilities to normalized ones (0-100%) @@ -609,7 +606,7 @@ def nn_predict_numpy(model, model_data, x): # apply input model normalisation x = (x - x_mean) / x_std - # gpu ? (usually not) + # gpu? if "current_gpu_device" not in model_data: current_gpu_mode, current_gpu_device = get_gpu_device(gpu_mode="auto") model_data["current_gpu_device"] = current_gpu_device @@ -716,7 +713,7 @@ def compute_angle_offset_torch(angles, length): cos_theta = torch.cos(angles_rad[:, 0]) cos_phi = torch.cos(angles_rad[:, 1]) - print('FIXME acos !! ') + print("FIXME acos !! ") # see in Gate_NN_ARF_Actor, line "phi = acos(dir.x())/degree;" tx = length * cos_phi @@ -914,14 +911,12 @@ def arf_from_points_to_image_counts( This version is model-agnostic (PyTorch or XGBoost). """ - # --- 1. Predict scatter probabilities with the correct model --- - # Get directions and energy for model input dirs = projected_batch[:, 2:5] energy = projected_batch[:, 5:6] # Calculate angles from directions - if model_data['angle_param'] == 'atan2': + if model_data["angle_param"] == "atan2": theta = np.degrees(np.arccos(np.clip(dirs[:, 2], -1, 1))) phi = np.degrees(np.arctan2(dirs[:, 1], dirs[:, 0])) else: @@ -945,8 +940,6 @@ def arf_from_points_to_image_counts( weights = projected_batch[:, 6] w_pred = w_pred * weights[:, np.newaxis] - # --- The rest of the function (geometric projection) remains the same --- - # Get initial position (px, py) cx = projected_batch[:, 0:2] From 644f551fdc31bfb461556fe017baa5ea8cb7624f Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Tue, 29 Jul 2025 16:36:48 +0200 Subject: [PATCH 18/20] Remove verbose flag from learning rate scheduler initialization. The `verbose` parameter was removed from `ReduceLROnPlateau` calls to simplify the scheduler configuration. This change ensures consistent and cleaner logging behavior. --- garf/garf_train.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/garf/garf_train.py b/garf/garf_train.py index f8336fe..5504754 100644 --- a/garf/garf_train.py +++ b/garf/garf_train.py @@ -74,7 +74,7 @@ def nn_get_optimiser(model_data, model): # decreasing learning_rate scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau( - optimizer, "min", verbose=False, patience=50 + optimizer, "min", patience=50 ) return optimizer, scheduler @@ -350,7 +350,7 @@ def train_nn(x_train, y_train, params): ) # decreasing learning_rate scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau( - optimizer, "min", verbose=False, patience=3 + optimizer, "min", patience=3 ) # Main loop initialization From 12b861e6ecceee80713d976f65d55eb3d368a404 Mon Sep 17 00:00:00 2001 From: David Sarrut Date: Mon, 25 Aug 2025 11:11:39 +0200 Subject: [PATCH 19/20] add debug bin --- garf/bin/garf_plot_training_dataset2.py | 89 ++++++++++++++++++++ garf/bin/garf_train_xgboost.py | 104 ++++++++++++++++++++++++ pyproject.toml | 1 + 3 files changed, 194 insertions(+) create mode 100755 garf/bin/garf_plot_training_dataset2.py create mode 100755 garf/bin/garf_train_xgboost.py diff --git a/garf/bin/garf_plot_training_dataset2.py b/garf/bin/garf_plot_training_dataset2.py new file mode 100755 index 0000000..c33454c --- /dev/null +++ b/garf/bin/garf_plot_training_dataset2.py @@ -0,0 +1,89 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import garf +import matplotlib.pyplot as plt +import numpy as np +import click + +# ----------------------------------------------------------------------------- +CONTEXT_SETTINGS = dict(help_option_names=["-h", "--help"]) + + +@click.command(context_settings=CONTEXT_SETTINGS) +@click.argument("data_file") +@click.option( + "--sample-size", default=5000, help="Number of non-detected (w=0) events to plot." +) +def garf_plot_training_dataset_2d(data_file, sample_size): + """ + \b + Display 2D scatter plots of the training dataset to show relationships + between angles, energy, and the final energy window. + + : dataset in root format + """ + print(f"Loading data from '{data_file}'") + data, theta, phi, E, w = garf.load_training_dataset(data_file) + print("Data loaded. Preparing plots...") + + # Find the unique window IDs and assign colors + window_ids = np.unique(w) + colors = plt.cm.viridis(np.linspace(0, 1, len(window_ids))) + + # Create a 1x3 figure for the plots + fig, ax = plt.subplots(1, 3, figsize=(18, 5.5)) + + # --- Sub-sample the dominant class for clarity --- + # Separate data by window ID + data_by_window = {win_id: data[w == win_id] for win_id in window_ids} + + # Sub-sample the window=0 class if it's too large + if 0 in data_by_window and len(data_by_window[0]) > sample_size: + print( + f"Sub-sampling non-detected (window=0) class from {len(data_by_window[0])} to {sample_size} points for clarity." + ) + indices = np.random.choice( + data_by_window[0].shape[0], sample_size, replace=False + ) + data_by_window[0] = data_by_window[0][indices] + + # --- Create the plots --- + plot_titles = ["Theta vs. Phi", "Energy vs. Theta", "Energy vs. Phi"] + plot_vars = [(1, 0), (0, 2), (1, 2)] # (x_col_idx, y_col_idx) from `data` array + + for i, p_ax in enumerate(ax): + for win_id, color in zip(window_ids, colors): + if win_id not in data_by_window: + continue + + subset = data_by_window[win_id] + x_data = subset[:, plot_vars[i][0]] + y_data = subset[:, plot_vars[i][1]] + + # For Energy plots, convert to keV + if plot_vars[i][1] == 2: # Energy is the Y-axis + y_data = y_data * 1000 + + p_ax.scatter( + x_data, + y_data, + color=color, + label=f"Window {int(win_id)}", + alpha=0.5, + s=5, + ) + + p_ax.set_title(plot_titles[i]) + p_ax.set_xlabel(f"{['Theta', 'Phi', 'Energy (keV)'][plot_vars[i][0]]}") + p_ax.set_ylabel(f"{['Theta', 'Phi', 'Energy (keV)'][plot_vars[i][1]]}") + p_ax.legend() + p_ax.grid(True, linestyle="--", alpha=0.6) + + plt.tight_layout() + plt.show() + + +# ----------------------------------------------------------------------------- +if __name__ == "__main__": + garf_plot_training_dataset_2d() diff --git a/garf/bin/garf_train_xgboost.py b/garf/bin/garf_train_xgboost.py new file mode 100755 index 0000000..85cd823 --- /dev/null +++ b/garf/bin/garf_train_xgboost.py @@ -0,0 +1,104 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import numpy as np +import click +import json +import joblib +import xgboost + +# Assuming 'garf' is a module in your project that contains 'load_training_dataset' +import garf + +# ----------------------------------------------------------------------------- +CONTEXT_SETTINGS = dict(help_option_names=["-h", "--help"]) + + +@click.command(context_settings=CONTEXT_SETTINGS) +@click.argument("json_params_file") +@click.argument("data_file") +@click.argument("output_model_file") +def garf_train_xgboost(json_params_file, data_file, output_model_file): + """ + \b + Train an XGBoost model for ARF based on a training dataset. + : Training parameters (H, L, rr, etc) in json format + : Dataset in root format + : Filename for the output trained model (e.g., 'model.joblib') + """ + + # --- 1. Load Data and Parameters --- + print(f"Loading data from '{data_file}'") + data, theta, phi, E, w = garf.load_training_dataset(data_file) + x_train = np.column_stack((theta, phi, E)) + y_train = w.astype(int) # Labels must be integers + # counts how many different unique windows? + n_ene_win = np.max(y_train) + 1 + print(f"Number of ENE win: {n_ene_win}") + + print(f"Loading parameters from '{json_params_file}'") + with open(json_params_file) as f: + params = json.load(f) + rr_factor = 500#params["RR"] FIXME + print(f"Russian Roulette factor found: {rr_factor}") + + # --- 2. Prepare Data and Metadata --- + print("Preparing data and calculating weights...") + + # Normalize input features (good practice, though XGBoost is less sensitive) + x_mean = np.mean(x_train, 0) + x_std = np.std(x_train, 0) + x_train_normalized = (x_train - x_mean) / x_std + + # Create model_data dictionary to save with the model + model_data = { + "x_mean": x_mean, + "x_std": x_std, + "rr": rr_factor, + "N": len(x_train), + "n_ene_win": n_ene_win, + "model_type": "xgboost", # Add model type for later use + } + + # Calculate sample weights using the robust hybrid method + # Note: XGBoost uses 'sample_weight' (one weight per training sample) + class_counts = np.bincount(y_train) + class_weights = 1.0 / (class_counts + 1e-9) + class_weights = np.ones_like(class_counts) + if rr_factor > 1 and len(class_weights) > 0: + class_weights[0] *= rr_factor + # class_weights = class_weights / np.mean(class_weights) + + # Create the final sample_weight array + sample_weights = class_weights[y_train] + print(f"Calculated class weights: {class_weights}") + + # --- 3. Define and Train XGBoost Model --- + print("Training XGBoost model...") + # These are starter parameters; you can tune them for better performance + model = xgboost.XGBClassifier( + objective="multi:softprob", # Output probabilities for each class + n_estimators=200, # Number of boosting rounds (trees) was 200 + max_depth=8, # Maximum depth of a tree + learning_rate=0.01, # Step size shrinkage + # use_label_encoder=False, + eval_metric="mlogloss", + n_jobs=-1, # Use all available CPU cores + ) + + # Train the model + model.fit(x_train_normalized, y_train, sample_weight=sample_weights) + print("Training complete.") + + # --- 4. Save the Model and Metadata --- + # We save both the trained model and the metadata needed for inference + output_data = {"model": model, "model_data": model_data} + + print(f"Saving trained model and metadata to '{output_model_file}'") + joblib.dump(output_data, output_model_file) + print("Done.") + + +# ----------------------------------------------------------------------------- +if __name__ == "__main__": + garf_train_xgboost() diff --git a/pyproject.toml b/pyproject.toml index 955e883..c06d393 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,6 +39,7 @@ dependencies = [ [project.scripts] garf_compare_image_profile = "garf.bin.garf_compare_image_profile:garf_compare_image_profile" garf_plot_training_dataset = "garf.bin.garf_plot_training_dataset:garf_plot_training_dataset" +garf_plot_training_dataset2 = "garf.bin.garf_plot_training_dataset:garf_plot_training_dataset2" garf_plot_test_dataset = "garf.bin.garf_plot_test_dataset:garf_plot_training_dataset" garf_train = "garf.bin.garf_train:garf_train" garf_nn_info = "garf.bin.garf_nn_info:garf_nn_info" From dfffee8aaaea7a38eab7a28ddaa630ffb288574a Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 25 Aug 2025 09:13:43 +0000 Subject: [PATCH 20/20] [pre-commit.ci] Automatic python formatting --- garf/bin/garf_train.py | 8 ++++++-- garf/bin/garf_train_xgboost.py | 2 +- garf/garf_train.py | 6 ++---- 3 files changed, 9 insertions(+), 7 deletions(-) diff --git a/garf/bin/garf_train.py b/garf/bin/garf_train.py index 00b85d4..7360e54 100755 --- a/garf/bin/garf_train.py +++ b/garf/bin/garf_train.py @@ -16,8 +16,12 @@ @click.argument("data") @click.argument("output") @click.option("--progress-bar/--no-progress-bar", default=True) -@click.option("--rr", default=None, help="RR value (overwrite the one in the param file)") -@click.option("--epoch", default=None, help="Nb of epoch (overwrite the one in the param file)") +@click.option( + "--rr", default=None, help="RR value (overwrite the one in the param file)" +) +@click.option( + "--epoch", default=None, help="Nb of epoch (overwrite the one in the param file)" +) def garf_train(param, data, output, rr, epoch, progress_bar): """ \b diff --git a/garf/bin/garf_train_xgboost.py b/garf/bin/garf_train_xgboost.py index 85cd823..82650c0 100755 --- a/garf/bin/garf_train_xgboost.py +++ b/garf/bin/garf_train_xgboost.py @@ -39,7 +39,7 @@ def garf_train_xgboost(json_params_file, data_file, output_model_file): print(f"Loading parameters from '{json_params_file}'") with open(json_params_file) as f: params = json.load(f) - rr_factor = 500#params["RR"] FIXME + rr_factor = 500 # params["RR"] FIXME print(f"Russian Roulette factor found: {rr_factor}") # --- 2. Prepare Data and Metadata --- diff --git a/garf/garf_train.py b/garf/garf_train.py index 5504754..59a9458 100644 --- a/garf/garf_train.py +++ b/garf/garf_train.py @@ -51,7 +51,7 @@ def nn_prepare_data(x_train, y_train, params): # this flag indicate that we use the new version of angle parametrisation # (acos + atan2 instead of acos + acos) - model_data["angle_param"] = 'atan2' + model_data["angle_param"] = "atan2" # copy param except comments for i in params: @@ -349,9 +349,7 @@ def train_nn(x_train, y_train, params): weight_decay=1e-4, ) # decreasing learning_rate - scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau( - optimizer, "min", patience=3 - ) + scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(optimizer, "min", patience=3) # Main loop initialization epoch_max = model_data["epoch_max"]