Merge pull request #74 from pim-n/dev

BUGFIX: incorrect shape comparison causing all saves to fail

CITATION.cff

@@ -0,0 +1,26 @@
+# This CITATION.cff file was generated with cffinit.
+# Visit https://bit.ly/cffinit to generate yours today!
+
+cff-version: 1.2.0
+title: PG-RAD - Primary Gamma RADiation Simulator
+message: >-
+  If you use this software, please cite it using the
+  metadata from this file.
+type: software
+authors:
+  - given-names: Pim
+    family-names: Nelissen
+    email: pi0274ne-s@student.lu.se
+    affiliation: Lund University
+repository-code: 'https://github.com/pim-n/pg-rad'
+abstract: >-
+  Primary Gamma RADiation (PG-RAD) is a deterministic
+  landscape simulator to simulate vehicle-based gamma source
+  localization scenarios.
+keywords:
+  - gamma spectrometry
+  - mobile gamma spectrometry
+  - emergency preparedness
+  - environmental radiation
+license: MIT
+

README.md

@@ -1,59 +1,32 @@
 [![Python 3.12](https://img.shields.io/badge/python-3.12-blue.svg)](https://www.python.org/downloads/release/python-312/)
 [![Tests](https://github.com/pim-n/pg-rad/actions/workflows/ci-tests.yml/badge.svg)](https://github.com/pim-n/pg-rad/actions/workflows/ci-tests.yml)
-# pg-rad
+# PG-RAD - Primary Gamma RADiation landscape simulator
-Primary Gamma RADiation landscape - Development
 
-## Clone
+PG-RAD is a command-line software package for simulating localisation of orphan sources using mobile gamma spectrometry. PG-RAD provides a framework for construction road geometries, arbitrary distributions of point sources, modelling the detector response to primary gammas. PG-RAD provides a configurable framework for constructing detector trajectories, defining radioactive source distributions, modelling detector response, and generating synthetic acquisition data. User input is specified through YAML configuration files, which makes simulations reproducible and easily shared.
-```
+
-git clone https://github.com/pim-n/pg-rad
+## Installation
-cd pg-rad
+
-git checkout dev
+PG-RAD is a Python 3 package and is only tested on x86_64 Linux systems. The following installation instructions were testing for a fresh virtual machine running Ubuntu 26.04 LTS.
+
+## Conda installation (Tested and recommended)
+
+1. Install git by `sudo apt update && sudo apt install git`
+2. Install miniforge by following [these](https://conda-forge.org/download/) instructions.
+3. Create a file called `environment.yml`, and paste the following in there:
+
+```yaml
+name: my-pgrad-env
+channels:
+  - conda-forge
+dependencies:
+  - python=3.12
+  - pip:
+    - git+ssh://git@github.com/pim-n/pg-rad.git@v0.1.0
 ```
+4. Run `conda env create -f environment.yml` to create the environment
+5. Run `conda activate pg-rad-analysis`. You should now be in the conda environment `my-pgrad-env`.
+6. To test if installation was succesful, run `pgrad --example --show`. If this runs without errors and produces visual output, PG-RAD is correctly installed.
 
-or
+## Manual installation
 
-```
+If you prefer another virtual environment, you still need git installed. Ensure the Python version of the environment is `>3.12.4` and `<3.13`. Then, with your virtual environment activated, run `pip install git+ssh://git@github.com/pim-n/pg-rad.git@main`. Run `pgrad --example --show` to check if the installation was successful.
-git@github.com:pim-n/pg-rad.git
-cd pg-rad
-git checkout dev
-```
-
-## Dependencies / venv
-
-With Python verion `>=3.12.4` and `<3.13`, create a virtual environment and install pg-rad.
-
-```
-python3 -m venv .venv
-source .venv/bin/activate
-```
-
-With the virtual environment activated, run:
-
-```
-pip install -e .[dev]
-```
-
-## Running example landscape
-
-The example landscape can be generated using the command-line interface. Still in the virtual environment, run
-
-```
-pgrad --test --loglevel DEBUG
-```
-
-## Tests
-
-Tests can be run with `pytest` from the root directory of the repository. With the virtual environment activated, run:
-
-```
-pytest
-```
-
-## Local viewing of documentation
-
-PG-RAD uses [Material for MkDocs](https://squidfunk.github.io/mkdocs-material/) for generating documentation. It can be locally viewed by (in the venv) running:
-```
-mkdocs serve
-```
-
-where you can add the `--livereload` flag to automatically update the documentation as you write to the Markdown files.

pyproject.toml

@@ -10,12 +10,12 @@ where = ["src"]
 "pg_rad.configs" = ["*.yml"]
 
 [project]
-name = "pg-rad"
+name = "pgrad"
-version = "0.2.1"
+version = "0.1.0"
 authors = [
   { name="Pim Nelissen", email="pi0274ne-s@student.lu.se" },
 ]
-description = "Primary Gamma RADiation Landscape"
+description = "PG-RAD - Primary Gamma RADiation Simulator"
 readme = "README.md"
 requires-python = ">=3.12.4,<3.13"
 dependencies = [
@@ -37,4 +37,4 @@ Issues = "https://github.com/pim-n/pg-rad/issues"
 
 [project.optional-dependencies]
 
-dev = ["pytest", "mkinit", "notebook", "mkdocs-material", "mkdocstrings-python", "mkdocs-jupyter", "flake8"]
+dev = ["pytest", "notebook", "mkdocs-material", "mkdocstrings-python", "mkdocs-jupyter", "flake8"]

src/pg_rad/background/background.py

@@ -11,15 +11,20 @@ def generate_background(
     cps_array: np.ndarray,
     detector: Detector,
     energy_keV: float,
+    lam_inp: int | None = None,
+    seed: int | None = None
 
 ) -> np.ndarray:
     """
     Generate synthetic background cps for a given detector and energy.
     """
+    if not lam_inp:
         ROI_lo, ROI_hi = get_roi_from_fwhm(detector, energy_keV)
         lam = get_cps_from_roi(detector, ROI_lo, ROI_hi)
+    else:
+        lam = lam_inp
 
-    rng = np.random.default_rng()
+    rng = np.random.default_rng(seed=seed)
     return rng.poisson(lam=lam, size=cps_array.shape)
 
 

src/pg_rad/data/angular_efficiencies/LU_HPGe_90.csv

@@ -18,6 +18,7 @@ angle,662,1173,1332
 160,1.113,1.096,1.099
 170,1.091,1.076,1.083
 180,1.076,1.066,1.078
+-180,1.076,1.066,1.078
 -170,1.102,1.091,1.093
 -160,1.122,1.100,1.102
 -150,1.128,1.105,1.093

src/pg_rad/data/angular_efficiencies/LU_NaIR.csv

@@ -0,0 +1,38 @@
+angle,662
+0,0.027
+10,0.162
+20,0.346
+30,0.517
+40,0.662
+50,0.794
+60,0.882
+70,0.947
+80,0.995
+90,1.000
+100,0.970
+110,0.895
+120,0.778
+130,0.649
+140,0.546
+150,0.477
+160,0.387
+170,0.267
+180,0.205
+-180,0.205
+-170,0.266
+-160,0.385
+-150,0.527
+-140,0.671
+-130,0.764
+-120,0.838
+-110,0.763
+-100,0.838
+-90,0.903
+-80,0.904
+-70,0.898
+-60,0.862
+-50,0.717
+-40,0.337
+-30,0.154
+-20,0.253
+-10,0.125

src/pg_rad/data/detectors.csv

@@ -2,3 +2,4 @@ name,type,is_isotropic
 dummy,NaI,true
 LU_NaI_3inch,NaI,true
 LU_HPGe_90,HPGe,false
+LU_NaIR,NaI,false

src/pg_rad/data/field_efficiencies/LU_NaIR.csv

@@ -0,0 +1,5 @@
+energy_keV,field_efficiency_m2
+0,0
+661.657,0.0261
+1173.228,0.0203
+1332.492,0.0166

src/pg_rad/inputparser/parser.py

@@ -104,7 +104,7 @@ class ConfigParser:
     def _parse_options(self) -> SimulationOptionsSpec:
         options = self.config.get("options", {})
 
-        allowed = {"air_density_kg_per_m3", "seed"}
+        allowed = {"air_density_kg_per_m3", "seed", "bkg_cps"}
         self._warn_unknown_keys(
             section="options",
             provided=set(options.keys()),
@@ -116,23 +116,33 @@ class ConfigParser:
             defaults.DEFAULT_AIR_DENSITY
         )
         seed = options.get("seed")
+        bkg_cps = options.get("bkg_cps")
 
         if not isinstance(air_density, float) or air_density <= 0:
             raise InvalidConfigValueError(
                 "options.air_density_kg_per_m3 must be a positive float "
                 "in kg/m^3."
             )
+
         if (
-            seed is not None or
+            seed is not None and
             (isinstance(seed, int) and seed <= 0)
         ):
             raise InvalidConfigValueError(
                 "Seed must be a positive integer value."
             )
 
+        if bkg_cps is not None and (
+            not isinstance(bkg_cps, int) or bkg_cps < 0
+        ):
+            raise InvalidConfigValueError(
+                "Background CPS must be an integer >= 0."
+            )
+
         return SimulationOptionsSpec(
             air_density=air_density,
             seed=seed,
+            bkg_cps=bkg_cps
         )
 
     def _parse_path(self) -> PathSpec:

src/pg_rad/inputparser/specs.py

@@ -18,6 +18,7 @@ class RuntimeSpec:
 class SimulationOptionsSpec:
     air_density: float
     seed: int | None = None
+    bkg_cps: int | None = None
 
 
 @dataclass

src/pg_rad/landscape/builder.py

@@ -107,7 +107,8 @@ class LandscapeBuilder:
 
     def set_point_sources(
         self,
-        *sources: AbsolutePointSourceSpec | RelativePointSourceSpec
+        *sources: AbsolutePointSourceSpec | RelativePointSourceSpec,
+        bounds_check: bool = False
     ):
         """Add one or more point sources to the world.
 
@@ -148,7 +149,7 @@ class LandscapeBuilder:
 
             # we dont support -x values, but negative y values are possible as
             # the path is centered in the y direction.
-            if not (
+            if bounds_check and not (
                 (0 <= pos[0] <= self._size[0]) and
                 (-0.5 * self._size[1] <= pos[1] <= 0.5 * self._size[1])
             ):

src/pg_rad/main.py

@@ -2,6 +2,7 @@ import argparse
 import logging
 import sys
 
+from numpy.random import SeedSequence
 from pandas.errors import ParserError
 
 from pg_rad.exceptions.exceptions import (
@@ -77,12 +78,22 @@ def main():
                 gamma_energy_keV: 661
 
         detector: LU_NaI_3inch
+
+        options:
+          seed: 1234
         """
     elif args.config:
         input_config = args.config
-
+    else:
+        logger.warning(
+            "No input provided. Try --example or --config path/to/config.yml. "
+        )
+        sys.exit(1)
     try:
         cp = ConfigParser(input_config).parse()
+        if cp.options.seed is None:
+            entr = SeedSequence().entropy
+            cp.options.seed = int(str(entr)[:6])
         landscape = LandscapeDirector.build_from_config(cp)
         output = SimulationEngine(
             landscape=landscape,
@@ -98,7 +109,7 @@ def main():
         if args.showplots:
             plotter.plot()
 
-        if not (args.save and args.showplots):
+        if not (args.save or args.showplots):
             logger.warning(
                 "No output produced. Use --save flag to save outputs and/or "
                 "--showplots to display interactive plots."

src/pg_rad/physics/fluence.py

@@ -108,7 +108,10 @@ def calculate_counts_along_path(
     landscape: "Landscape",
     detector: "Detector",
     velocity: float,
+    t_acq: float,
     points_per_segment: int = 10,
+    bkg_cps_input: int | None = None,
+    seed: int | None = None
 ) -> Tuple[np.ndarray, np.ndarray]:
     """Compute the counts recorded in each acquisition period in the landscape.
 
@@ -116,6 +119,7 @@ def calculate_counts_along_path(
         landscape (Landscape): _description_
         detector (Detector): _description_
         points_per_segment (int, optional): _description_. Defaults to 100.
+        bkg_cps_input (int | None, optional): Optional background CPS.
 
     Returns:
         Tuple[np.ndarray, np.ndarray]: Array of acquisition points and
@@ -148,17 +152,41 @@ def calculate_counts_along_path(
         landscape, full_positions, detector
     )
 
+    if bkg_cps_input is None:
         bkg = generate_background(
-        cps, detector, landscape.point_sources[0].isotope.E
+            cps, detector, landscape.point_sources[0].isotope.E,
+            seed=seed
+        )
+    elif bkg_cps_input == 0:
+        bkg = bkg_cps_input
+    else:
+        bkg = generate_background(
+            cps, detector, landscape.point_sources[0].isotope.E,
+            lam_inp=bkg_cps_input,
+            seed=seed
         )
 
     cps_with_bg = cps + bkg
-    # reshape so each segment is on a row
-    cps_per_seg = cps_with_bg.reshape(num_segments, points_per_segment)
 
-    du = s[1] - s[0]
+    # Integrate along time dimension. see e.g. Bukartas (2021 doccomp.)
-    integrated_counts = np.trapezoid(cps_per_seg, dx=du, axis=1) / velocity
+    t = s / abs(velocity)
-    int_counts_result = np.zeros(num_points)
-    int_counts_result[1:] = integrated_counts
 
-    return original_distances, s, cps_with_bg, int_counts_result, np.mean(bkg)
+    # acquisition bins
+    t_bins = np.arange(0, t[-1] + t_acq, t_acq)
+
+    int_counts = np.zeros(len(t_bins) - 1)
+    acq_points = np.zeros(len(t_bins) - 1)
+
+    for i in range(len(t_bins) - 1):
+        mask = (t >= t_bins[i]) & (t < t_bins[i + 1])
+
+        int_counts[i] = np.trapezoid(cps_with_bg[mask], t[mask])
+        acq_points[i] = np.mean(s[mask])
+
+    return (
+        acq_points[:-1],
+        s[:-1],
+        cps_with_bg[:-1],
+        int_counts[:-1],
+        np.mean(bkg)
+        )

src/pg_rad/plotting/landscape_plotter.py

@@ -7,6 +7,7 @@ from matplotlib.patches import Circle
 from pg_rad.landscape.landscape import Landscape
 
 logger = logging.getLogger(__name__)
+plt.set_loglevel(level='warning')
 
 
 class LandscapeSlicePlotter:

src/pg_rad/plotting/result_plotter.py

@@ -10,6 +10,9 @@ from pg_rad.simulator.outputs import SimulationOutput
 from pg_rad.landscape.landscape import Landscape
 
 
+plt.set_loglevel(level='warning')
+
+
 class ResultPlotter:
     def __init__(self, landscape: Landscape, output: SimulationOutput):
         self.landscape = landscape
@@ -98,13 +101,26 @@ class ResultPlotter:
     def _draw_counts(self, ax):
         x = self.count_rate_res.distance[1:]
         y = self.count_rate_res.integrated_counts[1:]
+
+        yerr = np.sqrt(y)
+        peak_idx = y.argmax()
+
+        ax.fill_between(
+            x,
+            y - yerr,
+            y + yerr,
+            color='red',
+            alpha=0.15,
+        )
+
         ax.plot(
-            x, y, color='r', linestyle='--',
+            x, y,
-            alpha=0.2, label=f'max(counts) = {y.max():.2f}'
+            color='r', linestyle='--', alpha=0.2,
+            label=rf'max(counts)={y[peak_idx]:.2f} $\pm$ {yerr[peak_idx]:.2f}'
         )
         ax.legend(handlelength=0, handletextpad=0, fancybox=True)
         ax.scatter(x, y, color='r', marker='x')
-        ax.set_title('Integrated counts')
+        ax.set_title('Counts')
         ax.set_xlabel('Arc length s [m]')
         ax.set_ylabel('N')
 
@@ -116,6 +132,7 @@ class ResultPlotter:
             ["Air density (kg/m^3)", round(self.landscape.air_density, 3)],
             ["Total path length (m)", round(self.landscape.path.length, 3)],
             ["Readout points", len(self.count_rate_res.integrated_counts)],
+            ["Seed", self.count_rate_res.seed],
             ["Mean background cps", round(self.count_rate_res.mean_bkg_cps, 3)]
         ]
 

src/pg_rad/simulator/engine.py

@@ -3,6 +3,7 @@ from typing import List
 from pg_rad.landscape.landscape import Landscape
 from pg_rad.simulator.outputs import (
     CountRateOutput,
+    DetectorOutput,
     SimulationOutput,
     SourceOutput
 )
@@ -31,10 +32,13 @@ class SimulationEngine:
 
         count_rate_results = self._calculate_count_rate_along_path()
         source_results = self._calculate_point_source_distance_to_path()
+        detector_results = self._generate_detector_output()
 
         return SimulationOutput(
             name=self.landscape.name,
+            size=self.landscape.size,
             count_rate=count_rate_results,
+            detector=detector_results,
             sources=source_results
         )
 
@@ -43,7 +47,10 @@ class SimulationEngine:
             calculate_counts_along_path(
                 self.landscape,
                 self.detector,
-                velocity=self.runtime_spec.speed
+                velocity=self.runtime_spec.speed,
+                bkg_cps_input=self.sim_spec.bkg_cps,
+                t_acq=self.runtime_spec.acquisition_time,
+                seed=self.sim_spec.seed
             )
         )
 
@@ -54,7 +61,8 @@ class SimulationEngine:
             sub_points,
             cps,
             int_counts,
-            mean_bkg_counts
+            mean_bkg_counts,
+            self.sim_spec.seed
         )
 
     def _calculate_point_source_distance_to_path(self) -> List[SourceOutput]:
@@ -79,3 +87,13 @@ class SimulationEngine:
             )
 
         return source_output
+
+    def _generate_detector_output(self) -> DetectorOutput:
+        return DetectorOutput(
+            name=self.detector.name,
+            type=self.detector.type,
+            is_isotropic=self.detector.is_isotropic,
+            field_eff=self.detector.get_efficiency(
+                self.landscape.point_sources[0].isotope.E
+            )
+        )

src/pg_rad/simulator/outputs.py

@@ -12,6 +12,7 @@ class CountRateOutput:
     cps: List[float]
     integrated_counts: List[float]
     mean_bkg_cps: List[float]
+    seed: int
 
 
 @dataclass
@@ -24,8 +25,18 @@ class SourceOutput:
     dist_from_path: float
 
 
+@dataclass
+class DetectorOutput:
+    name: str
+    type: str
+    is_isotropic: bool
+    field_eff: float
+
+
 @dataclass
 class SimulationOutput:
     name: str
+    size: tuple
+    detector: DetectorOutput
     count_rate: CountRateOutput
     sources: List[SourceOutput]

src/pg_rad/utils/export.py

@@ -1,15 +1,27 @@
+from dataclasses import asdict
 from datetime import datetime as dt
+import json
 import os
 import logging
 import re
 
-from numpy import array, full_like
+from numpy import array, full_like, ndarray, bool_
 from pandas import DataFrame
 
 from pg_rad.simulator.outputs import SimulationOutput
 
 
 logger = logging.getLogger(__name__)
+logging.getLogger("PIL").setLevel(logging.WARNING)
+
+
+class NumpyEncoder(json.JSONEncoder):
+    def default(self, obj):
+        if isinstance(obj, ndarray):
+            return obj.tolist()
+        elif isinstance(obj, bool_):
+            return bool(obj)
+        return super().default(obj)
 
 
 def generate_folder_name(sim: SimulationOutput) -> str:
@@ -32,12 +44,29 @@ def save_results(sim: SimulationOutput, folder_name: str) -> None:
         if ans.lower() == 'n':
             return
 
+    logger.debug(
+        f"Integrated counts: {list(sim.count_rate.integrated_counts)}"
+    )
+    logger.debug(
+        f"Distances: {list(sim.count_rate.distance)}"
+    )
+
     df = generate_df(sim)
     csv_name = generate_csv_name(sim)
     df.to_csv(f"{folder_name}/{csv_name}.csv", index=False)
+    param_dict = generate_sim_param_dict(sim)
+    with open(f"{folder_name}/parameters.json", 'w') as f:
+        json.dump(param_dict, f, cls=NumpyEncoder)
     logger.info(f"Simulation output saved to {folder_name}!")
 
 
+def generate_sim_param_dict(sim: SimulationOutput) -> dict:
+    """Parse simulation parameters and hyperparameters to dictionary."""
+    d = asdict(sim)
+    d.pop('count_rate')
+    return d
+
+
 def generate_df(sim: SimulationOutput) -> DataFrame:
     """Parse simulation output to CSV format and the name of CSV."""
 
@@ -46,10 +75,21 @@ def generate_df(sim: SimulationOutput) -> DataFrame:
         sim.count_rate.mean_bkg_cps
     )
 
+    east_coords = sim.count_rate.x[1:]
+    north_coords = sim.count_rate.y[1:]
+
+    if len(east_coords) != sim.count_rate.integrated_counts.shape[0]:
+        east_coords = None
+        north_coords = None
+        logger.warning(
+            "PG-RAD currently does not support this experimental path."
+            " Only ROI_P, ROI_BR and Dist will be saved."
+        )
+
     result_df = DataFrame(
         {
-            "East": sim.count_rate.x,
+            "East": east_coords,
-            "North": sim.count_rate.y,
+            "North": north_coords,
             "ROI_P": sim.count_rate.integrated_counts,
             "ROI_BR": br_array,
             "Dist": sim.count_rate.distance
@@ -62,13 +102,13 @@ def generate_df(sim: SimulationOutput) -> DataFrame:
 def generate_csv_name(sim: SimulationOutput) -> str:
     """Generate CSV name according to Alex' specification"""
     num_src = len(sim.sources)
+    src_ids = [str(i+1) for i in range(num_src)]
     bkg_cps = round(sim.count_rate.mean_bkg_cps)
     source_param_strings = [
         [
             str(round(s.activity))+"MBq",
             str(round(s.dist_from_path))+"m",
-            str(round(s.position[0])),
+            str(round(s.position[0]))+'_'+str(round(s.position[1]))
-            str(round(s.position[1])),
         ]
         for s in sim.sources
     ]
@@ -82,5 +122,6 @@ def generate_csv_name(sim: SimulationOutput) -> str:
 
         src_str = "_".join(src_str_array.flat)
 
-    csv_name = f"{num_src}_src_{bkg_cps}_cps_bkg_{src_str}"
+    src_ids_str = "_".join(src_ids)
+    csv_name = f"{src_ids_str}_src_{bkg_cps}_cps_bkg_{src_str}"
     return csv_name