remove mkinit requirement

prepare for release 0.1.0

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
+

pyproject.toml

@@ -10,12 +10,12 @@ where = ["src"]
 "pg_rad.configs" = ["*.yml"]
 
 [project]
-name = "pg-rad"
-version = "0.2.1"
+name = "pgrad"
+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.
     """
-    ROI_lo, ROI_hi = get_roi_from_fwhm(detector, energy_keV)
-    lam = get_cps_from_roi(detector, ROI_lo, ROI_hi)
+    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

@@ -1,4 +1,5 @@
 name,type,is_isotropic
 dummy,NaI,true
 LU_NaI_3inch,NaI,true
-LU_HPGe_90,HPGe,false
+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
     )
 
-    bkg = generate_background(
-        cps, detector, landscape.point_sources[0].isotope.E
-    )
+    if bkg_cps_input is None:
+        bkg = generate_background(
+            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]
-    integrated_counts = np.trapezoid(cps_per_seg, dx=du, axis=1) / velocity
-    int_counts_result = np.zeros(num_points)
-    int_counts_result[1:] = integrated_counts
+    # Integrate along time dimension. see e.g. Bukartas (2021 doccomp.)
+    t = s / abs(velocity)
 
-    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='--',
-            alpha=0.2, label=f'max(counts) = {y.max():.2f}'
+            x, y,
+            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,22 @@ 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:
+        east_coords = None
+        north_coords = None
+        logger.warning(
+            "PG-RAD currently does not support interpolation"
+            " of experimental paths for export. Only ROI_P, ROI_BR and Dist"
+            " will be saved."
+        )
+
     result_df = DataFrame(
         {
-            "East": sim.count_rate.x,
-            "North": sim.count_rate.y,
+            "East": east_coords,
+            "North": north_coords,
             "ROI_P": sim.count_rate.integrated_counts,
             "ROI_BR": br_array,
             "Dist": sim.count_rate.distance
@@ -62,6 +103,7 @@ 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 = [
         [
@@ -81,5 +123,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