Compare commits

..

6 Commits

27 changed files with 954 additions and 336 deletions

View File

@ -1,4 +0,0 @@
---
title: pg_rad.landscape.create_landscape_from_path
---
::: pg_rad.landscape.create_landscape_from_path

View File

@ -1,4 +0,0 @@
---
title: pg_rad.landscape.Landscape
---
::: pg_rad.landscape.Landscape

View File

@ -1,5 +0,0 @@
---
title: pg_rad.objects.Object
---
::: pg_rad.objects.Object

View File

@ -1,4 +0,0 @@
---
title: pg_rad.path.Path
---
::: pg_rad.path.Path

View File

@ -1,5 +0,0 @@
---
title: pg_rad.path.path_from_RT90
---
::: pg_rad.path.path_from_RT90

View File

@ -1,4 +0,0 @@
---
title: pg_rad.path.simplify_path
---
::: pg_rad.path.simplify_path

View File

@ -1,5 +0,0 @@
---
title: pg_rad.sources.PointSource
---
::: pg_rad.sources.PointSource

218
docs/config-spec.md Normal file
View File

@ -0,0 +1,218 @@
!!! note
To get started quickly, you may copy and modify one of the example configs found [here](quickstart.md#example-configs).
The config file must be a [YAML](https://yaml.org/) file. YAML is a serialization language that works with key-value pairs, but in a syntax more readable than some other alternatives. In YAML, the indentation matters. I
The remainder of this chapter will explain the different required and optionals keys, what they represent, and allowed values.
## Required keys
### Simulation options
The first step is to name the simulation, and define the speed of the vehicle (assumed constant) and acquisition time.
#### Landscape name
The name is a string, which may include spaces, numbers and special characters.
Examples:
```yaml
name: test_landscape
```
```yaml
name: Test Landscape 1
```
#### Acquisition time
The acquisition time of the detector in seconds.
Example:
```yaml
acquisition_time: 1
```
!!! note
All units in the config file must be specified in SI units, e.g. meters and seconds, unless the key contains a unit itself (e.g. `activity_MBq` means activity in MegaBequerels).
#### Vehicle speed
The speed of the vehicle in m/s. Currently, the vehicle speed must be assumed constant. An example could be
```yaml
speed: 13.89 # this is approximately 50 km/h
```
!!! note
The text after the `#` signifies a comment. PG-RAD will ignore this, but it can be helpful for yourself to write notes.
### Path
The `path` keyword is used to create a path for the detector to travel along. There are two ways to specify a path; from experimental data or by specifying a procedural path.
#### Path - Experimental data
Currently the only supported coordinate format is the RT90 (East, North) coordinate system. If you have experimental data in CSV format with columns for these coordinates, then you can load that path into PG-RAD as follows:
```yaml
path:
file: path/to/experimental_data.csv
east_col_name: East
north_col_name: North
```
#### Path - Procedural path
Alternatively, you can let PG-RAD generate a path for you. A procedural path can be specified with at least two subkeys: `length` and `segments`.
Currently supported segments are: `straight`, `turn_left` and `turn_right`, and are provided in a list under the `segments` subkey as follows:
```yaml
path:
segments:
- straight
- turn_left
- straight
```
The length must also be specified, using the `length` subkey. `length` can be specified in two ways: a list with the same length as the `segments` list
```yaml
path:
segments:
- straight
- turn_left
- straight
length:
- 500
- 250
- 500
```
which will assign that length (meters) to each segment. Alternatively, a single number can be passed:
```yaml
path:
segments:
- straight
- turn_left
- straight
length: 1250
```
Setting the length for the total path will cause PG-RAD to *randomly assign* portions of the total length to each segment.
Finally, there is also an option to specify the turn angle in degrees:
```yaml
path:
segments:
- straight
- turn_left: 90
- straight
length: 1250
```
Like with the lengths, if a turn segment has no angle specified, a random one (within pre-defined limits) will be taken.
!!! warning
Letting PG-RAD randomly assign lengths and angles can cause (expected) issues. That is because of physics restrictions. If the combination of length, angle (radius) and velocity of the vehicle is such that the centrifugal force makes it impossible to take this turn, PG-RAD will raise an error. To fix it, you can 1) reduce the speed; 2) define a smaller angle for the turn; or 3) assign more length to the turn segment.
!!! info
For more information about how procedural roads are generated, including the random sampling of lengths and angles, see X
### Sources
Currently, the only type of source supported is a point source. Point sources can be added under the `sources` key, where the **subkey is the name** of the source:
```yaml
sources:
my_source: ...
```
the source name should not contain spaces or special characters other than `_` or `-`. There are three required subkeys under `sources.my_source`, which are: `activity_MBq`, `isotope` and `position`.
#### Source activity
The source activity is in MegaBequerels and must be a strictly positive number:
```yaml
sources:
my_source:
activity_MBq: 100
```
#### Source isotope
The isotope for the point source. This must be a string, following the naming convention of the symbol followed by the number of nucleons, e.g. `Cs137`:
```yaml
sources:
my_source:
activity_MBq: 100
isotope: Cs137
```
!!! info
Currently the following isotopes are supported: `Cs137`
#### Source position
There are two ways to specify the source position. Either with absolute (x,y,z) coordinates
```yaml
sources:
my_source:
activity_MBq: 100
isotope: Cs137
position: [0, 0, 0]
```
or relative to the path, using the subkeys `along_path`, `dist_from_path` and `side`
```yaml
sources:
my_source:
activity_MBq: 100
isotope: Cs137
position:
along_path: 100
dist_from_path: 50
side: left
```
Note that side is relative to the direction of travel. The path will by default start at (x,y) = (0,0) and initial heading is parallel to the x-axis.
### Detector
The final required key is the `detector`. Currently, only isotropic detectors are supported. Nonetheless, you must specify it with `name`, `is_isotropic` and `efficiency`:
```yaml
detector:
name: test
is_isotropic: True
efficiency: 0.02
```
Note there are some existing detectors available, where efficiency is not required and will be looked up by PG-RAD itself:
```yaml
detector:
name: NaIR
is_isotropic: True
```
## Optional keys
The following subkeys are optional and should be put under the `options` key.
```yaml
options:
air_density_kg_per_m3: 1.243
seed: 1234
```

File diff suppressed because one or more lines are too long

View File

@ -0,0 +1,118 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "1a063d05",
"metadata": {},
"source": [
"# Pseudo-ramdom procedural roads\n",
"\n",
"Suppose one wishes to describe a road between A and B in terms of segments\n",
"\n",
"$$\n",
"\\text{straight, turn left, straight, turn right, straight}\n",
"$$\n",
"\n",
"Let's see how we can create a random road of length $L$ from a pre-determined set of prefabs.\n",
"\n",
"#### Random apportionment of total length\n",
"\n",
"Suppose we want to build a road of length $L$ out of $K$ segments. The total number of waypoints $N$ depends on the step size $\\Delta s$:\n",
"\n",
"$$\n",
"N = \\frac{L}{\\Delta s}.\n",
"$$\n",
"\n",
"Let $\\left( p_1, p_2, \\dots, p_K \\right)$ represent the proportion of $N$ that each prefab will be assigned, where $\\sum p_i = 1$. One useful distribution here is the [Dirichlet distribution](https://en.wikipedia.org/wiki/Dirichlet_distribution), which is parametrized by a vector $\\mathbf{\\alpha} = \\left(\\alpha_1, \\alpha_2, \\dots, \\alpha_K \\right)$. The special case where all $\\alpha_i$, the scalar parameter $\\alpha$ is called a *concentration parameter*. Setting the same $\\alpha$ across the entire parameter space makes the distribution symmetric, meaning no prior assumptions are made regarding the proportion of $N$ that will be assigned to each segment. $\\alpha = 1$ leads to what is known as a flat Dirichlet distribution, whereas higher values lead to more dense and evenly distributed $\\left( p_1, p_2, \\dots, p_K \\right)$. On the other hand, keeping $\\alpha \\leq 1$ gives a sparser distribution which can lead to larger variance in apportioned number of waypoints to $\\left( p_1, p_2, \\dots, p_K \\right)$.\n",
"\n",
"#### Expectation value and variance of Dirichlet distribution\n",
"\n",
"Suppose we draw our samples for proportion of length from the Dirichlet distribution\n",
"\n",
"$$\n",
"(p_1, p_2, \\ldots, p_K) \\sim \\text{Dirichlet}(\\alpha, \\alpha, \\ldots, \\alpha)\n",
"$$\n",
"\n",
"with $\\alpha _{0}=\\sum _{i=1}^{K}\\alpha _{i}$, the mean and variance are then\n",
"\n",
"$$\n",
"\\operatorname {E} [p_{i}]={\\frac {\\alpha _{i}}{\\alpha _{0}}}, \\; \\operatorname {Var} [p_{i}]={\\frac {\\alpha _{i}(\\alpha _{0}-\\alpha _{i})}{\\alpha _{0}^{2}(\\alpha _{0}+1)}}.\n",
"$$\n",
"\n",
"If $\\alpha$ is a scalar, then $\\alpha _{0}= K \\alpha$ and the above simplifies to\n",
"\n",
"$$\n",
"\\operatorname {E} [p_{i}]={\\frac {\\alpha}{K \\alpha}}={\\frac {1}{K}}, \\; \\operatorname {Var} [p_{i}]={\\frac {\\alpha(K \\alpha -\\alpha)}{(K \\alpha)^{2}(K \\alpha +1)}}.\n",
"$$\n",
"\n",
"We see that $\\operatorname {Var} [p_{i}] \\propto \\frac{1}{\\alpha}$ meaning that the variance reduces with increasing $\\alpha$. We can simply scale the proportions\n",
"\n",
"$$\n",
"(N \\cdot p_1, N \\cdot p_2, \\ldots, N \\cdot p_K)\n",
"$$\n",
"\n",
"to get the randomly assigned number of waypoints for each prefab. We now have a distribution which can give randomly assigned lengths to a given list of prefabs, with a parameter to control the degree of randomness. With a large concentration parameter $\\alpha$, the distribution of lengths will be more uniform, with each prefab getting $N \\cdot \\operatorname {E} [p_{i}]={\\frac {N}{K}}$ waypoints assigned to it. Likewise, keeping $\\alpha$ low increases variance and allows for a more random assignment of proportions of waypoints to each prefab segment.\n",
"\n",
"#### Random angles\n",
"\n",
"Suppose a turn of a pre-defined arc length $l$ made of $N/K$ waypoints. If one wants to create a random angle, one has to keep in mind that the minimum radius $R_{min}$ depends on the speed of the vehicle and the weather conditions:\n",
"\n",
"$$\n",
"R_{\\text{min,vehicle}} = \\frac{v^2}{g\\mu},\n",
"$$\n",
"\n",
"where\n",
"- $v$ is the velocity of the vehicle in $\\text{m/s}$,\n",
"- $g$ is the gravitational acceleration (about $9.8$ $\\text{m/s}^{2}$), and\n",
"- $\\mu$ is the friction coefficient (about $0.7$ for dry asphalt).\n",
"\n",
"A regular turn (not a U-turn or roundabout) should also have an lower and upper limit on the angle, say, 30 degrees to 90 degrees for a conservative estimate. In terms of radii, it becomes\n",
"\n",
"$$\n",
"R_{\\text{min}} = \\max\\left(R_{\\text{min,vehicle}}, \\frac{l}{\\pi/2}\\right)\n",
"$$\n",
"\n",
"and\n",
"\n",
"$$\n",
"R_{\\text{max}} = \\frac{l}{\\pi/6}.\n",
"$$\n",
"\n",
"We then sample\n",
"\n",
"$$\n",
"R \\sim \\text{Uniform}\\left(R_{\\text{min}}, R_{\\text{max\\_angle}}\\right)\n",
"$$\n",
"\n",
"and obtain a random radius for a turn of arc length $l$ with limits to ensure the radius is large enough given the velocity of the vehicle. Finally, the curvature profile is related to the radius by\n",
"\n",
"$$\n",
"\\kappa = \\frac{1}{R}\n",
"$$\n",
"\n",
"which means that the curvature profile of a turn is simply a vector $\\mathbf{\\kappa} = (1/R, \\dots, 1/R)$ with a length of $N/K$ waypoints."
]
}
],
"metadata": {
"kernelspec": {
"display_name": ".venv",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.9"
}
},
"nbformat": 4,
"nbformat_minor": 5
}

View File

@ -2,29 +2,17 @@
Primary Gamma RADiation Landscapes (PG-RAD) is a Python package for research in source localization. It can simulate mobile gamma spectrometry data acquired from vehicle-borne detectors along a predefined path (e.g. a road). Primary Gamma RADiation Landscapes (PG-RAD) is a Python package for research in source localization. It can simulate mobile gamma spectrometry data acquired from vehicle-borne detectors along a predefined path (e.g. a road).
## Requirements ## About
PG-RAD requires Python `3.12`. The guides below assume a unix-like system. This software has been developed as part of dissertation work for the degree of master of Computational Science and Physics at Lund University, Sweden. The work has been done at the department of Medical Radiation Physics (MSF), Faculty of Medicine. The radiological emergency preparedness research group of MSF is assigned by the Swedish Radiation Safety Authority (SSM) to aid in preparation for effective mitigation of radiological or nuclear disasters on Swedish soil.
## Installation (CLI) ## Value proposition
<!--pipx seems like a possible option to install python package in a contained environment on unix--> PG-RAD is a toolbox that allows for simulation of detector response for a wide variety of source localization scenarios. The strength of the software lies in its simple and minimal configuration and user input, while its flexibility allows for reconstruction of specific scenarios with relative ease. PG-RAD is also general enough that novel methods such as UAV-borne detectors can be simulated and evaluated.
Lorem ipsum User input takes the form of an input file (YAML), describing the path, detector and source(s), and optional parameters. The output of the program is visualizations of the world (the path and sources), as well as the detector count rate as a function of distance travelled along the path.
## Installation (Python module) Users can provide experimental / geographical coordinates representing real roads. Alternatively, users can let PG-RAD generate a procedural road, where the user can easily control what that road should look like. The user can specify a single point source, several point sources, as well as a field of radioactive material covering a large area.
If you are interested in using PG-RAD in another Python project, create a virtual environment first:
```
python3 -m venv .venv
```
Then install PG-RAD in it:
```
source .venv/bin/activate
(.venv) pip install git+https://github.com/pim-n/pg-rad
``` ```
See how to get started with PG-RAD with your own Python code [here](pg-rad-in-python). See how to get started with PG-RAD with your own Python code [here](pg-rad-in-python).

47
docs/installation.md Normal file
View File

@ -0,0 +1,47 @@
## Requirements
PG-RAD requires Python `>=3.12.4` and `<3.13`. It has been tested on `3.12.9`. The guides below assume a unix-like system. You can check the Python version you have installed as follows:
```
python --version
```
If you don't have the right version installed there are various ways to get a compatible version, such as [pyenv](https://github.com/pyenv/pyenv?tab=readme-ov-file#installation).
## Installation (CLI)
<!--pipx seems like a possible option to install python package in a contained environment on unix-->
Lorem ipsum
## Installation (Python module)
If you are interested in using PG-RAD in another Python project, create a virtual environment first:
```
python -m venv .venv
```
Then install PG-RAD in it:
```
source .venv/bin/activate
(.venv) pip install git+https://github.com/pim-n/pg-rad
```
See how to get started with PG-RAD with your own Python code [here](pg-rad-in-python).
## For developers
```
git clone https://github.com/pim-n/pg-rad
cd pg-rad
git checkout dev
```
or
```
git@github.com:pim-n/pg-rad.git
cd pg-rad
git checkout dev
```

View File

@ -1,4 +0,0 @@
---
title: Using PG-RAD in CLI
---
Lorem ipsum.

File diff suppressed because one or more lines are too long

187
docs/quickstart.md Normal file
View File

@ -0,0 +1,187 @@
## Installation
See the [installation guide](installation.md).
## Test your installation
First, see if PG-RAD is available on your system by typing
```
pgrad --help
```
You should get output along the lines of
```
usage: pg-rad [-h] ...
Primary Gamma RADiation landscape tool
...
```
If you get something like `pgrad: command not found`, please consult the [installation guide](installation.md).
You can run a quick test scenario as follows:
```
pgrad --test
```
This should produce a plot of a scenario containing a single point source and a path.
## Running PG-RAD
In order to use the CLI for your own simulations, you need to provide a *config file*. To run with your config, run
```
pgrad --config path/to/my_config.yml
```
where `path/to/my_config.yml` points to your config file.
## Example configs
The easiest way is to take one of these example configs, and adjust them as needed. Alternatively, there is a detailed guide on how to write your own config file [here](config-spec.md).
=== "Example 1"
The position can be defined relative to the path. `along_path` means at what distance traveled along the path the source is found. If the path is 200 meters long and `along_path` is `100` then the source is halfway along the path. `dist_from_path` is the distance in meters from the path. `side` is the side of the path the source is located. This is relative to the direction the path is traveled.
``` yaml
name: Example 1
speed: 13.89
acquisition_time: 1
path:
file: path/to/exp_coords.csv
east_col_name: East
north_col_name: North
sources:
source1:
activity_MBq: 1000
isotope: CS137
position:
along_path: 100
dist_from_path: 50
side: left
detector:
name: dummy
is_isotropic: True
```
=== "Example 2"
The position can also just be defined with (x,y,z) coordinates.
``` yaml
name: Example 2
speed: 13.89
acquisition_time: 1
path:
file: path/to/exp_coords.csv
east_col_name: East
north_col_name: North
sources:
source1:
activity_MBq: 1000
isotope: CS137
position: [104.3, 32.5, 0]
source2:
activity_MBq: 100
isotope: CS137
position: [0, 0, 0]
detector:
name: dummy
is_isotropic: True
```
=== "Example 3"
This is an example of a procedural path with random apportionment of total length and random angles being assigned to turns. The parameter `alpha` is optional, and is related to randomness. A higher value leads to more uniform apportionment of lengths and a lower value to more random apportionment. More information about `alpha` can be found [here](pg-rad-config-spec.md).
``` yaml
name: Example 3
speed: 8.33
acquisition_time: 1
path:
length: 1000
segments:
- straight
- turn_left
- straight
alpha: 100
sources:
source1:
activity_MBq: 1000
isotope: CS137
position: [0, 0, 0]
detector:
name: dummy
is_isotropic: True
```
=== "Example 4"
This is an example of a procedural path that is partially specified. Note that turn_left now is a key for the corresponding angle of 45 degrees. The length is still divided randomly
``` yaml
name: Example 4
speed: 8.33
acquisition_time: 1
path:
length: 1000
segments:
- straight
- turn_left: 45
- straight
sources:
source1:
activity_MBq: 1000
isotope: CS137
position: [0, 0, 0]
detector:
name: dummy
is_isotropic: True
```
=== "Example 5"
This is an example of a procedural path that is fully specified. See how length is now a list matching the length of the segments.
``` yaml
name: Example 5
speed: 8.33
acquisition_time: 1
path:
length:
- 400
- 200
- 400
segments:
- straight
- turn_left: 45
- straight
sources:
source1:
activity_MBq: 1000
isotope: CS137
position: [0, 0, 0]
detector:
name: dummy
is_isotropic: True
```

View File

@ -30,6 +30,11 @@ markdown_extensions:
- pymdownx.superfences - pymdownx.superfences
- pymdownx.arithmatex: - pymdownx.arithmatex:
generic: true generic: true
- admonition
- pymdownx.details
- pymdownx.tabbed:
alternate_style: true
combine_header_slug: true
extra_javascript: extra_javascript:
- javascripts/mathjax.js - javascripts/mathjax.js
@ -47,3 +52,12 @@ plugins:
options: options:
show_source: false show_source: false
show_root_heading: false show_root_heading: false
nav:
- Home: index.md
- Installation Guide: installation.md
- Quickstart Guide: quickstart.md
- 'Tutorial: Writing a Config File': config-spec.md
- Explainers:
- explainers/planar_curve.ipynb
- explainers/prefab_roads.ipynb

View File

@ -3,10 +3,13 @@ disable_existing_loggers: false
formatters: formatters:
simple: simple:
format: '%(asctime)s - %(levelname)s: %(message)s' format: '%(asctime)s - %(levelname)s: %(message)s'
colored:
'()': pg_rad.logger.logger.ColorFormatter
format: '%(asctime)s - %(levelname)s: %(message)s'
handlers: handlers:
stdout: stdout:
class: logging.StreamHandler class: logging.StreamHandler
formatter: simple formatter: colored
stream: ext://sys.stdout stream: ext://sys.stdout
loggers: loggers:
root: root:

View File

@ -14,7 +14,7 @@ class DetectorBuilder:
if self.detector_spec.is_isotropic: if self.detector_spec.is_isotropic:
return IsotropicDetector( return IsotropicDetector(
self.detector_spec.name, self.detector_spec.name,
self.detector_spec.eff self.detector_spec.efficiency
) )
else: else:
raise NotImplementedError("Angular detector not supported yet.") raise NotImplementedError("Angular detector not supported yet.")

View File

@ -5,10 +5,10 @@ class BaseDetector(ABC):
def __init__( def __init__(
self, self,
name: str, name: str,
eff: float efficiency: float
): ):
self.name = name self.name = name
self.eff = eff self.efficiency = efficiency
def get_efficiency(self): def get_efficiency(self):
pass pass
@ -18,21 +18,21 @@ class IsotropicDetector(BaseDetector):
def __init__( def __init__(
self, self,
name: str, name: str,
eff: float | None = None efficiency: float,
): ):
super().__init__(name, eff) super().__init__(name, efficiency)
def get_efficiency(self, energy): def get_efficiency(self, energy):
return self.eff return self.efficiency
class AngularDetector(BaseDetector): class AngularDetector(BaseDetector):
def __init__( def __init__(
self, self,
name: str, name: str,
eff: float | None = None efficiency: float
): ):
super().__init__(name, eff) super().__init__(name, efficiency)
def get_efficiency(self, angle, energy): def get_efficiency(self, angle, energy):
pass pass

View File

@ -31,3 +31,7 @@ class DimensionError(ValueError):
class InvalidIsotopeError(ValueError): class InvalidIsotopeError(ValueError):
"""Raised if attempting to load an isotope that is not valid.""" """Raised if attempting to load an isotope that is not valid."""
class InvalidConfigValueError(ValueError):
"""Raised if a config key has an incorrect type or value."""

View File

@ -4,7 +4,11 @@ from typing import Any, Dict, List, Union
import yaml import yaml
from pg_rad.exceptions.exceptions import MissingConfigKeyError, DimensionError from pg_rad.exceptions.exceptions import (
MissingConfigKeyError,
DimensionError,
InvalidConfigValueError
)
from pg_rad.configs import defaults from pg_rad.configs import defaults
from .specs import ( from .specs import (
@ -14,7 +18,7 @@ from .specs import (
PathSpec, PathSpec,
ProceduralPathSpec, ProceduralPathSpec,
CSVPathSpec, CSVPathSpec,
SourceSpec, PointSourceSpec,
AbsolutePointSourceSpec, AbsolutePointSourceSpec,
RelativePointSourceSpec, RelativePointSourceSpec,
DetectorSpec, DetectorSpec,
@ -98,19 +102,23 @@ class ConfigParser:
def _parse_options(self) -> SimulationOptionsSpec: def _parse_options(self) -> SimulationOptionsSpec:
options = self.config.get("options", {}) options = self.config.get("options", {})
allowed = {"air_density", "seed"} allowed = {"air_density_kg_per_m3", "seed"}
self._warn_unknown_keys( self._warn_unknown_keys(
section="options", section="options",
provided=set(options.keys()), provided=set(options.keys()),
allowed=allowed, allowed=allowed,
) )
air_density = options.get("air_density", defaults.DEFAULT_AIR_DENSITY) air_density = options.get(
"air_density_kg_per_m3",
defaults.DEFAULT_AIR_DENSITY
)
seed = options.get("seed") seed = options.get("seed")
if not isinstance(air_density, float) or air_density <= 0: if not isinstance(air_density, float) or air_density <= 0:
raise ValueError( raise ValueError(
"options.air_density must be a positive float in kg/m^3." "options.air_density_kg_per_m3 must be a positive float "
"in kg/m^3."
) )
if ( if (
seed is not None or seed is not None or
@ -230,9 +238,9 @@ class ConfigParser:
def _is_turn(segment_type: str) -> bool: def _is_turn(segment_type: str) -> bool:
return segment_type in {"turn_left", "turn_right"} return segment_type in {"turn_left", "turn_right"}
def _parse_point_sources(self) -> List[SourceSpec]: def _parse_point_sources(self) -> List[PointSourceSpec]:
source_dict = self.config.get("sources", {}) source_dict = self.config.get("sources", {})
specs: List[SourceSpec] = [] specs: List[PointSourceSpec] = []
for name, params in source_dict.items(): for name, params in source_dict.items():
@ -245,7 +253,7 @@ class ConfigParser:
isotope = params.get("isotope") isotope = params.get("isotope")
if not isinstance(activity, int | float) or activity <= 0: if not isinstance(activity, int | float) or activity <= 0:
raise ValueError( raise InvalidConfigValueError(
f"sources.{name}.activity_MBq must be positive value " f"sources.{name}.activity_MBq must be positive value "
"in MegaBequerels." "in MegaBequerels."
) )
@ -270,6 +278,16 @@ class ConfigParser:
) )
elif isinstance(position, dict): elif isinstance(position, dict):
alignment = position.get("acquisition_alignment")
if alignment not in {'best', 'worst', None}:
raise InvalidConfigValueError(
f"sources.{name}.acquisition_alignment must be "
"'best' or 'worst', with 'best' aligning source "
f"{name} in the middle of the two nearest acquisition "
"points, and 'worst' aligning exactly perpendicular "
"to the nearest acquisition point."
)
specs.append( specs.append(
RelativePointSourceSpec( RelativePointSourceSpec(
name=name, name=name,
@ -278,12 +296,13 @@ class ConfigParser:
along_path=float(position["along_path"]), along_path=float(position["along_path"]),
dist_from_path=float(position["dist_from_path"]), dist_from_path=float(position["dist_from_path"]),
side=position["side"], side=position["side"],
z=position.get("z", defaults.DEFAULT_SOURCE_HEIGHT) z=position.get("z", defaults.DEFAULT_SOURCE_HEIGHT),
alignment=alignment
) )
) )
else: else:
raise ValueError( raise InvalidConfigValueError(
f"Invalid position format for source '{name}'." f"Invalid position format for source '{name}'."
) )
@ -295,7 +314,7 @@ class ConfigParser:
missing = required - det_dict.keys() missing = required - det_dict.keys()
if missing: if missing:
raise MissingConfigKeyError(missing) raise MissingConfigKeyError("detector", missing)
name = det_dict.get("name") name = det_dict.get("name")
is_isotropic = det_dict.get("is_isotropic") is_isotropic = det_dict.get("is_isotropic")
@ -303,19 +322,23 @@ class ConfigParser:
default_detectors = defaults.DETECTOR_EFFICIENCIES default_detectors = defaults.DETECTOR_EFFICIENCIES
if eff is None and name in default_detectors.keys(): if name in default_detectors.keys() and not eff:
eff = default_detectors[name] eff = default_detectors[name]
elif eff is not None: elif eff:
pass pass
else: else:
raise ValueError( raise ValueError(
f"The detector {name} is not in the library, and no " f"The detector {name} not found in library. Either "
"efficiency was defined. Either specify detector efficiency " f"specify {name}.efficiency or "
"or choose one from the following list: " "choose a detector from the following list: "
f"{default_detectors.keys()}" f"{default_detectors.keys()}."
) )
return DetectorSpec(name=name, eff=eff, is_isotropic=is_isotropic) return DetectorSpec(
name=name,
efficiency=eff,
is_isotropic=is_isotropic
)
def _warn_unknown_keys(self, section: str, provided: set, allowed: set): def _warn_unknown_keys(self, section: str, provided: set, allowed: set):
unknown = provided - allowed unknown = provided - allowed

View File

@ -1,5 +1,6 @@
from abc import ABC from abc import ABC
from dataclasses import dataclass from dataclasses import dataclass
from typing import Literal
@dataclass @dataclass
@ -42,31 +43,32 @@ class CSVPathSpec(PathSpec):
@dataclass @dataclass
class SourceSpec(ABC): class PointSourceSpec(ABC):
activity_MBq: float activity_MBq: float
isotope: str isotope: str
name: str name: str
@dataclass @dataclass
class AbsolutePointSourceSpec(SourceSpec): class AbsolutePointSourceSpec(PointSourceSpec):
x: float x: float
y: float y: float
z: float z: float
@dataclass @dataclass
class RelativePointSourceSpec(SourceSpec): class RelativePointSourceSpec(PointSourceSpec):
along_path: float along_path: float
dist_from_path: float dist_from_path: float
side: str side: str
z: float z: float
alignment: Literal["best", "worst"] | None
@dataclass @dataclass
class DetectorSpec: class DetectorSpec:
name: str name: str
eff: float | None efficiency: float
is_isotropic: bool is_isotropic: bool
@ -76,5 +78,5 @@ class SimulationSpec:
runtime: RuntimeSpec runtime: RuntimeSpec
options: SimulationOptionsSpec options: SimulationOptionsSpec
path: PathSpec path: PathSpec
point_sources: list[SourceSpec] point_sources: list[PointSourceSpec]
detector: DetectorSpec detector: DetectorSpec

View File

@ -40,7 +40,7 @@ class CS137(Isotope):
preset_isotopes: Dict[str, Type[Isotope]] = { preset_isotopes: Dict[str, Type[Isotope]] = {
"CS137": CS137 "Cs137": CS137
} }

View File

@ -1,5 +1,7 @@
import logging import logging
from typing import Self from typing import Literal, Self
import numpy as np
from .landscape import Landscape from .landscape import Landscape
from pg_rad.dataloader.dataloader import load_data from pg_rad.dataloader.dataloader import load_data
@ -113,11 +115,25 @@ class LandscapeBuilder:
pos = (s.x, s.y, s.z) pos = (s.x, s.y, s.z)
elif isinstance(s, RelativePointSourceSpec): elif isinstance(s, RelativePointSourceSpec):
path = self.get_path() path = self.get_path()
if s.alignment:
along_path = self._align_relative_source(
s.along_path,
path,
s.alignment
)
logger.info(
f"Because source {s.name} was set to align with path "
f"({s.alignment} alignment), it was moved to be at "
f"{along_path} m along the path from {s.along_path} m."
)
else:
along_path = s.along_path
pos = rel_to_abs_source_position( pos = rel_to_abs_source_position(
x_list=path.x_list, x_list=path.x_list,
y_list=path.y_list, y_list=path.y_list,
path_z=path.z, path_z=path.z,
along_path=s.along_path, along_path=along_path,
side=s.side, side=s.side,
dist_from_path=s.dist_from_path) dist_from_path=s.dist_from_path)
if any( if any(
@ -158,6 +174,48 @@ class LandscapeBuilder:
max_size = max(self._path.size) max_size = max(self._path.size)
self.set_landscape_size((max_size, max_size)) self.set_landscape_size((max_size, max_size))
def _align_relative_source(
self,
along_path: float,
path: "Path",
mode: Literal["best", "worst"],
) -> tuple[float, float, float]:
"""Given the arc length at which the point source is placed,
align the source relative to the waypoints of the path. Here,
'best' means the point source is moved such that it is
perpendicular to the midpoint between two acuisition points.
'worst' means the point source is moved such that it is
perpendicular to the nearest acquisition point.
The distance to the path is not affected by this algorithm.
For more details on alignment, see
Fig. 4 and page 24 in Bukartas (2021).
Args:
along_path (float): Current arc length position of the source.
path (Path): The path to align to.
mode (Literal["best", "worst"]): Alignment mode.
Returns:
along_new (float): The updated arc length position.
"""
ds = np.hypot(
path.x_list[1] - path.x_list[0],
path.y_list[1] - path.y_list[0],
)
if mode == "worst":
along_new = round(along_path / ds) * ds
elif mode == "best":
along_new = (round(along_path / ds - 0.5) + 0.5) * ds
else:
raise ValueError(f"Unknown alignment mode: {mode}")
return along_new
def build(self): def build(self):
landscape = Landscape( landscape = Landscape(
name=self.name, name=self.name,

View File

@ -22,7 +22,7 @@ class LandscapeDirector:
def build_test_landscape(): def build_test_landscape():
fp = files('pg_rad.data').joinpath(TEST_EXP_DATA) fp = files('pg_rad.data').joinpath(TEST_EXP_DATA)
source = PointSource( source = PointSource(
activity_MBq=100E9, activity_MBq=100E6,
isotope="CS137", isotope="CS137",
position=(0, 0, 0) position=(0, 0, 0)
) )

View File

@ -20,3 +20,21 @@ def setup_logger(log_level: str = "WARNING"):
config["loggers"]["root"]["level"] = log_level config["loggers"]["root"]["level"] = log_level
logging.config.dictConfig(config) logging.config.dictConfig(config)
class ColorFormatter(logging.Formatter):
# ANSI escape codes
COLORS = {
logging.DEBUG: "\033[36m", # Cyan
logging.INFO: "\033[32m", # Green
logging.WARNING: "\033[33m", # Yellow
logging.ERROR: "\033[31m", # Red
logging.CRITICAL: "\033[41m", # Red background
}
RESET = "\033[0m"
def format(self, record):
color = self.COLORS.get(record.levelno, self.RESET)
record.levelname = f"{color}{record.levelname}{self.RESET}"
record.msg = f"{record.msg}"
return super().format(record)

View File

@ -10,6 +10,7 @@ from pg_rad.exceptions.exceptions import (
MissingConfigKeyError, MissingConfigKeyError,
OutOfBoundsError, OutOfBoundsError,
DimensionError, DimensionError,
InvalidConfigValueError,
InvalidIsotopeError InvalidIsotopeError
) )
from pg_rad.logger.logger import setup_logger from pg_rad.logger.logger import setup_logger
@ -64,15 +65,20 @@ def main():
activity_MBq: 1000 activity_MBq: 1000
position: [500, 100, 0] position: [500, 100, 0]
isotope: CS137 isotope: CS137
detector:
name: dummy
is_isotropic: True
""" """
cp = ConfigParser(test_yaml).parse() cp = ConfigParser(test_yaml).parse()
landscape = LandscapeDirector.build_from_config(cp) landscape = LandscapeDirector.build_from_config(cp)
detector = DetectorBuilder(cp.detector).build()
output = SimulationEngine( output = SimulationEngine(
landscape=landscape, landscape=landscape,
runtime_spec=cp.runtime, runtime_spec=cp.runtime,
sim_spec=cp.options sim_spec=cp.options,
detector=detector
).simulate() ).simulate()
plotter = ResultPlotter(landscape, output) plotter = ResultPlotter(landscape, output)
@ -96,18 +102,19 @@ def main():
MissingConfigKeyError, MissingConfigKeyError,
KeyError, KeyError,
YAMLError, YAMLError,
): ) as e:
logger.critical(e)
logger.critical( logger.critical(
"The provided config file is invalid. " "The config file is missing required keys or may be an "
"Check the log above. You can consult the documentation for " "invalid YAML file. Check the log above. Consult the "
"an explanation of how to define a config file." "documentation for examples of how to write a config file."
) )
sys.exit(1) sys.exit(1)
except ( except (
OutOfBoundsError, OutOfBoundsError,
DimensionError, DimensionError,
InvalidIsotopeError, InvalidIsotopeError,
ValueError InvalidConfigValueError
) as e: ) as e:
logger.critical(e) logger.critical(e)
logger.critical( logger.critical(