From 01461965d9daad30cd296d5594f230157fca4b70 Mon Sep 17 00:00:00 2001 From: Niklas Wahl Date: Wed, 5 Aug 2026 10:46:31 +0200 Subject: [PATCH 1/6] Add documentation --- .github/workflows/docs.yml | 54 ++++ .gitignore | 5 + .readthedocs.yaml | 16 ++ CMakeLists.txt | 19 ++ docs/Doxyfile | 60 +++++ docs/README.md | 23 ++ docs/c-api/engines.md | 21 ++ docs/c-api/geometry-and-sources.md | 19 ++ docs/c-api/index.md | 19 ++ docs/c-api/infrastructure.md | 25 ++ docs/conf.py | 185 +++++++++++++ docs/getting-started/installation.md | 29 ++ docs/index.md | 56 ++++ docs/python-api/index.md | 76 ++++++ docs/requirements.txt | 12 + pyproject.toml | 11 + src/omc_engine_cube.h | 93 ++++--- src/omc_engine_dij.h | 68 +++-- src/omc_engine_forward.h | 106 ++++---- src/omc_geom.h | 61 +++-- src/omc_host.h | 79 ++++-- src/omc_random.h | 76 ++++-- src/omc_score.h | 121 +++++---- src/omc_source_beamlet.h | 78 +++--- src/omc_spectrum.h | 56 ++-- src/omc_utilities.h | 195 +++++++++----- ucodes/omc_python/ompmc/__init__.py | 387 +++++++++++++++++++++++---- 27 files changed, 1531 insertions(+), 419 deletions(-) create mode 100644 .github/workflows/docs.yml create mode 100644 .readthedocs.yaml create mode 100644 docs/Doxyfile create mode 100644 docs/README.md create mode 100644 docs/c-api/engines.md create mode 100644 docs/c-api/geometry-and-sources.md create mode 100644 docs/c-api/index.md create mode 100644 docs/c-api/infrastructure.md create mode 100644 docs/conf.py create mode 100644 docs/getting-started/installation.md create mode 100644 docs/index.md create mode 100644 docs/python-api/index.md create mode 100644 docs/requirements.txt diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml new file mode 100644 index 0000000..c2a8069 --- /dev/null +++ b/.github/workflows/docs.yml @@ -0,0 +1,54 @@ +name: docs + +# Builds the Sphinx/Doxygen/Breathe site the same way Read the Docs does (see +# .readthedocs.yaml), so a broken build is caught on the PR rather than after +# it merges. Deliberately does not install the project itself -- the docs +# build needs no C++ compiler, see the comment in docs/conf.py -- only +# docs/requirements.txt and the doxygen binary. + +on: + push: + branches: ['**'] + pull_request: + paths: + - 'docs/**' + - 'src/**' + - 'ucodes/omc_python/**' + - 'README.md' + - 'BUILDING.md' + - 'CMakeLists.txt' + - '.readthedocs.yaml' + - '.github/workflows/docs.yml' + workflow_dispatch: + +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + +jobs: + docs: + runs-on: ubuntu-24.04 + steps: + - uses: actions/checkout@v4 + + - uses: actions/setup-python@v5 + with: + python-version: '3.12' + + - name: Install doxygen + run: sudo apt-get update && sudo apt-get install -y doxygen + + - name: Install docs dependencies + run: pip install -r docs/requirements.txt + + # -W turns warnings into errors, -n reports broken cross-references + # (a :func:`...` to a symbol that no longer exists, for example), and + # --keep-going collects every failure in one run instead of stopping at + # the first. + - name: Build docs + run: sphinx-build -b html -W --keep-going -n docs docs/_build/html + + - uses: actions/upload-artifact@v4 + with: + name: docs-html + path: docs/_build/html diff --git a/.gitignore b/.gitignore index 99ca69e..c52ee8b 100644 --- a/.gitignore +++ b/.gitignore @@ -18,6 +18,11 @@ __pycache__/ dist/ .pytest_cache/ +# Sphinx/Doxygen documentation build +docs/_build/ +docs/doxygen/ +.venv-docs/ + # CMake build trees and their products build/ build-*/ diff --git a/.readthedocs.yaml b/.readthedocs.yaml new file mode 100644 index 0000000..80f0048 --- /dev/null +++ b/.readthedocs.yaml @@ -0,0 +1,16 @@ +version: 2 + +build: + os: ubuntu-24.04 + tools: + python: "3.12" + apt_packages: + - doxygen + commands: + # `build.commands` replaces RTD's default build entirely, so this does + # its own dependency install (docs/requirements.txt, NOT the project -- + # see the comment there) and its own sphinx-build. -W/-n turn broken + # cross-references and undocumented autodoc targets into a failed build, + # matching .github/workflows/docs.yml. + - pip install -r docs/requirements.txt + - sphinx-build -b html -W --keep-going -n docs $READTHEDOCS_OUTPUT/html diff --git a/CMakeLists.txt b/CMakeLists.txt index 82e4f95..9f9eaee 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -733,6 +733,25 @@ if(OMPMC_BUILD_TESTS) endif() endif() +################################################################################ +# Documentation +################################################################################ + +# Convenience target for people who drive everything through CMake; not part +# of ALL and not needed for Read the Docs, which runs sphinx-build directly +# (see .readthedocs.yaml) so that building the docs never needs this +# configure step or a C++ toolchain. See docs/README.md for the equivalent +# plain sphinx-build command. +find_program(SPHINX_BUILD_EXECUTABLE sphinx-build) +if(SPHINX_BUILD_EXECUTABLE) + add_custom_target(docs + COMMAND "${SPHINX_BUILD_EXECUTABLE}" -b html + "${CMAKE_CURRENT_SOURCE_DIR}/docs" + "${CMAKE_CURRENT_BINARY_DIR}/docs/html" + WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}" + COMMENT "Building Sphinx documentation") +endif() + ################################################################################ # Installation ################################################################################ diff --git a/docs/Doxyfile b/docs/Doxyfile new file mode 100644 index 0000000..edb7a91 --- /dev/null +++ b/docs/Doxyfile @@ -0,0 +1,60 @@ +# Minimal Doxygen config: XML output only, consumed by Breathe/Sphinx. Runs as +# a subprocess of docs/conf.py, not part of the CMake build -- see the comment +# there for why. Any setting not listed here uses Doxygen's built-in default. + +PROJECT_NAME = "ompMC" +PROJECT_NUMBER = $(OMPMC_DOXYGEN_VERSION) +PROJECT_BRIEF = "OpenMP parallel Monte Carlo photon and electron transport" +OUTPUT_DIRECTORY = doxygen +QUIET = YES +WARN_AS_ERROR = NO + +# Public C API only: the 10 headers a host embeds ompMC through. ompmc.h/.c +# hold the physics transport internals, meant for maintainers reading the +# source directly rather than an API reference; see docs/c-api/internals.md. +INPUT = ../src +FILE_PATTERNS = *.h +EXCLUDE = ../src/ompmc.h +RECURSIVE = NO + +# The GPL banner opens with /**** rather than /** or /*!, so Doxygen skips it +# without an EXCLUDE_PATTERNS entry -- see JAVADOC_AUTOBRIEF below. +JAVADOC_AUTOBRIEF = YES +OPTIMIZE_OUTPUT_FOR_C = YES +TYPEDEF_HIDES_STRUCT = NO +EXTRACT_ALL = NO +EXTRACT_STATIC = NO +HIDE_UNDOC_MEMBERS = NO +HIDE_UNDOC_CLASSES = NO +SORT_MEMBER_DOCS = NO +MARKDOWN_SUPPORT = YES +WARN_IF_UNDOCUMENTED = YES + +# Each header's @file comment mentions the others by plain filename (e.g. +# "(omc_geom.h)") as prose, not as a link. With autolinking on, Doxygen turns +# those into cross-file references that only resolve if every referenced +# file's page happens to already be registered in the same Sphinx build in +# the right order, which our per-file `doxygenfile` pages spread across three +# separate c-api/*.md pages do not reliably satisfy -- breaking the build +# with "undefined label" errors that have nothing to do with a real missing +# doc. Explicit references (`#name`, `struct Foo`) are unaffected; this only +# turns off the automatic kind. +AUTOLINK_SUPPORT = NO + +# omcLog()/omcFail() carry GCC/MSVC attribute macros (noreturn, printf format +# checking) that Doxygen's C parser copies verbatim into the declaration text +# it hands Breathe, and Breathe's declaration parser then chokes on them ("not +# valid C++") since they are not real syntax. Expanding just these two to +# nothing, rather than turning on full preprocessing, keeps every other macro +# in the headers (MXMED, MAX_MEDIA-style constants, threadprivate guards) +# showing up in the docs as written. +ENABLE_PREPROCESSING = YES +MACRO_EXPANSION = YES +EXPAND_ONLY_PREDEF = YES +PREDEFINED = OMC_NORETURN= \ + OMC_PRINTF_LIKE(fmtArg,firstArg)= + +GENERATE_HTML = NO +GENERATE_LATEX = NO +GENERATE_XML = YES +XML_PROGRAMLISTING = NO diff --git a/docs/README.md b/docs/README.md new file mode 100644 index 0000000..41f1b92 --- /dev/null +++ b/docs/README.md @@ -0,0 +1,23 @@ +# Building the documentation locally + +No C++ compiler needed: the compiled `_ompmc` extension is mocked out (see +`conf.py`), so this only needs Python, Doxygen, and the packages in +`requirements.txt`. + +```sh +# Doxygen: apt install doxygen / brew install doxygen / choco install doxygen.strawberry +python -m venv .venv-docs +.venv-docs/Scripts/activate # .venv-docs/bin/activate on Linux/macOS +pip install -r docs/requirements.txt + +sphinx-build -b html docs docs/_build/html +``` + +Open `docs/_build/html/index.html`. Rebuild after an edit with the same +command; add `-E` to force a full rebuild (e.g. after editing `conf.py` or +`Doxyfile`) or `-W -n --keep-going` to build the way CI does, which turns +warnings and broken cross-references into a failing build. + +`sphinx-autobuild docs docs/_build/html` (from the `sphinx-autobuild` +package, not in `requirements.txt`) rebuilds and reloads a browser tab on +every save, useful while editing prose. diff --git a/docs/c-api/engines.md b/docs/c-api/engines.md new file mode 100644 index 0000000..6ef61db --- /dev/null +++ b/docs/c-api/engines.md @@ -0,0 +1,21 @@ +# Dose engines + +The three ways to run a calculation, differing only in where particles start +and how the result comes back. All three share the phantom +({doc}`geometry-and-sources`), the spectrum and the transport itself, and are +singletons: one calculation at a time per process. + +## Dose-influence matrix (Dij) + +```{doxygenfile} omc_engine_dij.h +``` + +## Forward dose from weighted beamlets + +```{doxygenfile} omc_engine_forward.h +``` + +## Dose cube from a collimated beam + +```{doxygenfile} omc_engine_cube.h +``` diff --git a/docs/c-api/geometry-and-sources.md b/docs/c-api/geometry-and-sources.md new file mode 100644 index 0000000..40e3db7 --- /dev/null +++ b/docs/c-api/geometry-and-sources.md @@ -0,0 +1,19 @@ +# Geometry and sources + +The rectilinear voxel phantom every engine transports in, and the two things +that turn it plus a spectrum into starting particles. + +## Voxel phantom + +```{doxygenfile} omc_geom.h +``` + +## Beamlet apertures + +```{doxygenfile} omc_source_beamlet.h +``` + +## Source spectrum + +```{doxygenfile} omc_spectrum.h +``` diff --git a/docs/c-api/index.md b/docs/c-api/index.md new file mode 100644 index 0000000..ae52855 --- /dev/null +++ b/docs/c-api/index.md @@ -0,0 +1,19 @@ +# C API + +`ompmc_core` is a C99 static library shared by every user code -- the +`omc_dosxyz` command line binary, the `omc_matrad` MATLAB/Octave MEX file, and +the `_ompmc` Python extension. A new host reaches it through ten headers, +grouped below the way {doc}`../getting-started/installation` introduces them. + +The transport physics itself (`ompmc.c`/`ompmc.h`, roughly 6,000 lines of +photon and electron interaction code) is not part of this reference: it has no +public entry points a host calls directly, and is meant to be read as source +rather than browsed as an API. + +```{toctree} +:maxdepth: 2 + +engines +geometry-and-sources +infrastructure +``` diff --git a/docs/c-api/infrastructure.md b/docs/c-api/infrastructure.md new file mode 100644 index 0000000..57022dc --- /dev/null +++ b/docs/c-api/infrastructure.md @@ -0,0 +1,25 @@ +# Infrastructure + +Cross-cutting pieces every user code and every engine shares: how shared code +talks back to its host, dose/variance scoring, random numbers, and small +utilities. + +## Host callbacks + +```{doxygenfile} omc_host.h +``` + +## Scoring + +```{doxygenfile} omc_score.h +``` + +## Random numbers + +```{doxygenfile} omc_random.h +``` + +## Utilities + +```{doxygenfile} omc_utilities.h +``` diff --git a/docs/conf.py b/docs/conf.py new file mode 100644 index 0000000..3174a52 --- /dev/null +++ b/docs/conf.py @@ -0,0 +1,185 @@ +"""Sphinx configuration for the ompMC documentation. + +Builds with no compiler and no CMake configure step: the C API comes from +Doxygen (invoked below as a subprocess) and the Python API from the pure +Python ``ompmc`` package under ``ucodes/omc_python``, imported straight from +source with its compiled ``_ompmc`` extension mocked out. See +docs/README.md for how to build this locally. +""" + +from __future__ import annotations + +import os +import re +import subprocess +import sys +from pathlib import Path + +DOCS_DIR = Path(__file__).resolve().parent +ROOT_DIR = DOCS_DIR.parent + +# So `import ompmc` finds the pure Python package without installing the +# compiled extension (see autodoc_mock_imports below). +sys.path.insert(0, str(ROOT_DIR / "ucodes" / "omc_python")) + +# -- Project information ----------------------------------------------------- + +project = "ompMC" +copyright = "2018, Edgardo Doerner" +author = "Edgardo Doerner" + + +def _version_from_cmake() -> str: + """The single source of truth for the version is CMakeLists.txt's + project() call; pyproject.toml reads it with the same regex.""" + text = (ROOT_DIR / "CMakeLists.txt").read_text(encoding="utf-8") + match = re.search(r"project\(ompMC\s+VERSION\s+([0-9]+\.[0-9]+\.[0-9]+)", text) + return match.group(1) if match else "0.0.0" + + +release = _version_from_cmake() +version = ".".join(release.split(".")[:2]) + +# -- General configuration ---------------------------------------------------- + +extensions = [ + "myst_parser", + "sphinx.ext.autodoc", + "sphinx.ext.napoleon", + "sphinx.ext.intersphinx", + "sphinx.ext.viewcode", + "sphinx.ext.autosectionlabel", + "breathe", + "sphinx_copybutton", + "sphinx_design", +] + +templates_path = [] +# README.md is about building the *docs*, for someone reading the repo +# source; it is not a page of the published site. +exclude_patterns = ["_build", "Thumbs.db", ".DS_Store", "doxygen", "README.md"] + +# BUILDING.md is included verbatim into getting-started/installation.md (see +# there) so it stays a single source of truth, but its GitHub-relative links +# to files outside docs/ (src/, LICENSE, .github/...) are not part of this +# site and MyST cannot resolve them as cross-references. That is a known, +# accepted gap -- see installation.md's note -- rather than a broken link in +# hand-authored content, which is why it is silenced globally instead of +# fixed link by link. +suppress_warnings = [ + "myst.xref_missing", + # Breathe's `doxygenfile` directive emits a self-referencing permalink + # for the file compound itself (distinct from the functions/structs + # inside it, which link fine); that target is never actually registered, + # so nitpicky mode (-n) reports it as broken on every c-api page + # regardless of what the file documents. Tracked upstream in breathe; + # nothing on our side to fix. + "ref.ref", +] + +source_suffix = { + ".rst": "restructuredtext", + ".md": "markdown", +} + +autosectionlabel_prefix_document = True + +# -- MyST ---------------------------------------------------------------- + +myst_enable_extensions = [ + "colon_fence", + "deflist", + "fieldlist", +] +myst_heading_anchors = 3 + +# -- Autodoc / Napoleon ---------------------------------------------------- + +# The compiled extension needs a C++ toolchain, CMake and OpenMP to build. +# Nothing in the public API requires it to be importable: the pure Python +# wrapper in ompmc/__init__.py is what is documented, so the compiled +# submodule is mocked instead of built. +autodoc_mock_imports = ["ompmc._ompmc"] + +autodoc_member_order = "bysource" +autodoc_typehints = "description" +autodoc_default_options = { + "members": True, + "undoc-members": False, + "show-inheritance": True, +} + +napoleon_google_docstring = False +napoleon_numpy_docstring = True +napoleon_use_rtype = False +napoleon_attr_annotations = True + +intersphinx_mapping = { + "python": ("https://docs.python.org/3", None), + "numpy": ("https://numpy.org/doc/stable/", None), + "scipy": ("https://docs.scipy.org/doc/scipy/", None), +} + +# -- Doxygen / Breathe --------------------------------------------------- + +DOXYGEN_XML_DIR = DOCS_DIR / "doxygen" / "xml" + + +def _run_doxygen() -> None: + """Regenerate the Doxygen XML that breathe reads from. + + Run unconditionally on every build (local or Read the Docs) rather than + from CMake, so that `sphinx-build docs docs/_build/html` is the entire + command and docs never depend on a configured CMake tree. Read the Docs + installs the `doxygen` binary via `build.apt_packages` in + .readthedocs.yaml. + """ + doxyfile = DOCS_DIR / "Doxyfile" + env = dict(os.environ, OMPMC_DOXYGEN_VERSION=release) + subprocess.run(["doxygen", str(doxyfile)], cwd=DOCS_DIR, env=env, check=True) + + +if os.environ.get("SKIP_DOXYGEN") != "1": + _run_doxygen() + +breathe_projects = {"ompMC": str(DOXYGEN_XML_DIR)} +breathe_default_project = "ompMC" +breathe_default_members = ("members", "undoc-members") +# Without this, breathe renders .h files under the C++ domain (a header could +# be either), and OPTIMIZE_OUTPUT_FOR_C in the Doxyfile is wasted. +breathe_domain_by_extension = {"h": "c"} + +# The C domain tries to cross-reference every type name it sees in a +# signature. 's fixed-width types are never going to be Doxygen +# output (they are not part of this project's INPUT), so nitpick mode (-n) +# would otherwise fail the build over them. +nitpick_ignore_regex = [ + ("c:identifier", r"u?int(8|16|32|64)_t"), + ("c:identifier", r"size_t"), + # Documented as c:macro entries (see omc_utilities.h), which the + # automatic xref generated for an array-size expression in a struct + # member's signature does not look up under. + ("c:identifier", r"BUFFER_SIZE"), + ("c:identifier", r"INPUT_PAIRS"), + # Napoleon renders a numpydoc "type, optional" field through the Python + # domain's TypedField, which cross-references every comma/"or"-separated + # token in it -- including words that describe the type rather than name + # one. These are the vocabulary used across ompmc/__init__.py's + # docstrings, not real classes. + ("py:class", r"optional"), + ("py:class", r"callable"), + ("py:class", r"array_like"), + ("py:class", r"sequence"), +] + +# -- HTML output ------------------------------------------------------------ + +html_theme = "furo" +html_static_path = ["_static"] +html_title = f"ompMC {version}" + +html_theme_options = { + "source_repository": "https://github.com/e0404/ompMC", + "source_branch": "master", + "source_directory": "docs/", +} diff --git a/docs/getting-started/installation.md b/docs/getting-started/installation.md new file mode 100644 index 0000000..6d01c74 --- /dev/null +++ b/docs/getting-started/installation.md @@ -0,0 +1,29 @@ +# Installation + +The Python package installs the ordinary way: + +```sh +pip install ompmc +``` + +which compiles the extension for the interpreter it is run with; a C++ +compiler and a working OpenMP runtime are all it needs. Prebuilt wheels for +Linux, macOS and Windows carry their own OpenMP runtime, so they need +neither. + +Building the C library, `omc_dosxyz` and the MATLAB/Octave MEX file uses +CMake directly. What follows is the build guide from the repository root, +included here so it has one copy. A few of its relative links (to `src/`, +`LICENSE`, and similar files outside `docs/`) resolve on GitHub but not from +this site; the [repository on GitHub](https://github.com/e0404/ompMC) is the +canonical place to browse those. +so it has one copy. A few of its relative links (to `src/`, `LICENSE`, and +similar files outside `docs/`) resolve on GitHub but not from this site; the +[repository on GitHub](https://github.com/e0404/ompMC) is the canonical place +to browse those. + +```{include} ../../BUILDING.md +:relative-docs: docs/ +:relative-images: +:start-line: 1 +``` diff --git a/docs/index.md b/docs/index.md new file mode 100644 index 0000000..bdf8c01 --- /dev/null +++ b/docs/index.md @@ -0,0 +1,56 @@ +# ompMC + +ompMC is an OpenMP-parallelized, CPU-based Monte Carlo code for coupled +photon-electron transport in voxelized geometries, built for beamlet-based +treatment planning: `omc_matrad` transports histories for many beamlets in one +run and returns the dose-influence matrix (Dij) an optimizer needs for +fluence-map optimization, and the Python interface gives the same +calculations to any other host. + +::::{grid} 1 2 2 2 +:gutter: 3 + +:::{grid-item-card} {octicon}`rocket` Getting started +:link: getting-started/installation +:link-type: doc +Build ompMC and its Python extension. +::: + +:::{grid-item-card} {octicon}`package` Python API +:link: python-api/index +:link-type: doc +`Geometry`, `Spectrum`, `calc_dij`, `calc_cube`, `calc_forward`. +::: + +:::{grid-item-card} {octicon}`file-code` C API +:link: c-api/index +:link-type: doc +The engine headers a new host embeds ompMC through. +::: + +:::{grid-item-card} {octicon}`mark-github` Source +:link: https://github.com/e0404/ompMC +The repository, issue tracker and full README. +::: +:::: + +```{toctree} +:hidden: +:caption: Getting started + +getting-started/installation +``` + +```{toctree} +:hidden: +:caption: Python + +python-api/index +``` + +```{toctree} +:hidden: +:caption: C + +c-api/index +``` diff --git a/docs/python-api/index.md b/docs/python-api/index.md new file mode 100644 index 0000000..526f2bb --- /dev/null +++ b/docs/python-api/index.md @@ -0,0 +1,76 @@ +# Python API + +The `ompmc` package wraps the compiled `_ompmc` extension: dataclasses +describe the phantom, source and physics, and three functions run a +calculation against them, sharing the same phantom, physics and source +spectra. + +```{list-table} +:header-rows: 1 + +* - Function + - Result +* - {py:func}`ompmc.calc_dij` + - One sparse column of dose per beamlet -- the dose-influence matrix a + treatment planning system optimizes against. +* - {py:func}`ompmc.calc_forward` + - Dose everywhere in the phantom from a whole weighted set of beamlets, in + one go -- what `calc_dij(...) @ weights` would give. +* - {py:func}`ompmc.calc_cube` + - Dose everywhere in the phantom from a single collimated beam. +``` + +## Quickstart + +```{literalinclude} ../../README.md +:language: python +:start-after: "```python" +:end-before: "```" +:dedent: +``` + +```{note} +Cubes must be Fortran ordered (see {py:class}`~ompmc.Geometry` below): the +transport indexes voxels with the first axis varying fastest, so a C-ordered +cube would describe a transposed phantom -- it is rejected rather than +silently copied. +``` + +## Phantom, sources and physics + +```{eval-rst} +.. autoclass:: ompmc.Geometry + :members: + +.. autoclass:: ompmc.Spectrum + :members: + +.. autoclass:: ompmc.BeamletSource + :members: + +.. autoclass:: ompmc.CollimatedSource + :members: + +.. autoclass:: ompmc.Physics + :members: +``` + +## Calculations + +```{eval-rst} +.. autofunction:: ompmc.calc_dij + +.. autofunction:: ompmc.calc_forward + +.. autofunction:: ompmc.calc_cube +``` + +## Utilities + +```{eval-rst} +.. autofunction:: ompmc.data_path + +.. data:: ompmc.MAX_MEDIA + + Maximum number of distinct media a :class:`Geometry` may reference. +``` diff --git a/docs/requirements.txt b/docs/requirements.txt new file mode 100644 index 0000000..3b70d65 --- /dev/null +++ b/docs/requirements.txt @@ -0,0 +1,12 @@ +# Documentation build dependencies. Deliberately NOT the project itself +# (which needs a C++ compiler, CMake and OpenMP to build its compiled +# extension) -- conf.py adds ucodes/omc_python to sys.path and mocks the +# compiled _ompmc submodule instead. numpy is a real import at module load +# time in ompmc/__init__.py, so it is listed here rather than mocked. +sphinx>=7 +furo +myst-parser +breathe +sphinx-copybutton +sphinx-design +numpy>=1.22 diff --git a/pyproject.toml b/pyproject.toml index f3e47d2..119e05b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,6 +24,17 @@ Homepage = "https://github.com/e0404/ompMC" [project.optional-dependencies] test = ["pytest >=7"] +# Listed for reference; the docs build itself uses docs/requirements.txt so +# that building the site never needs the compiled extension -- see the +# comment in docs/conf.py. +docs = [ + "sphinx >=7", + "furo", + "myst-parser", + "breathe", + "sphinx-copybutton", + "sphinx-design", +] [tool.scikit-build] # 3.26 is where FindPython grew the Development.SABIModule component that the diff --git a/src/omc_engine_cube.h b/src/omc_engine_cube.h index 71041a5..e893a4b 100644 --- a/src/omc_engine_cube.h +++ b/src/omc_engine_cube.h @@ -19,7 +19,8 @@ along with this program. If not, see . *****************************************************************************/ -/****************************************************************************** +/*! + @file omc_engine_cube - Dose in every voxel from one collimated beam. This is what omc_dosxyz calculates: a point source at a given distance from @@ -38,8 +39,10 @@ 3. built the source spectrum (omc_spectrum.h) 4. called omcSsdSourceInit() on the source below - and afterwards it owns the cleanup of those. Like the rest of ompMC this is a - singleton: one calculation at a time per process. + and afterwards it owns the cleanup of those. + + @warning Like the rest of ompMC this is a singleton: one calculation at a + time per process. *****************************************************************************/ #ifndef OMC_ENGINE_CUBE_H @@ -49,69 +52,85 @@ struct OmcSpectrum; -/* A point source at distance ssd in front of the phantom, shining through a +/*! A point source at distance ssd in front of the phantom, shining through a rectangular collimator opening on the phantom surface. The host fills the first five members; omcSsdSourceInit() clamps the rectangle to the phantom and works out the rest. */ struct OmcSsdSource { - double ssd; // distance of point source to phantom surface + double ssd; ///< distance of point source to phantom surface - double xinl, xinu; // lower and upper x-bounds of the field on - // phantom surface - double yinl, yinu; // lower and upper y-bounds of the field on - // phantom surface + double xinl; ///< lower x-bound of the field on the phantom surface + double xinu; ///< upper x-bound of the field on the phantom surface + double yinl; ///< lower y-bound of the field on the phantom surface + double yinu; ///< upper y-bound of the field on the phantom surface /* Derived by omcSsdSourceInit() */ - double xsize, ysize; // x- and y-width of collimated field - int ixinl, ixinu; // lower and upper x-bounds indices of the - // field on phantom surface - int iyinl, iyinu; // lower and upper y-bounds indices of the - // field on phantom surface + double xsize; ///< x-width of the collimated field + double ysize; ///< y-width of the collimated field + int ixinl; ///< voxel index of the field's lower x-bound + int ixinu; ///< voxel index of the field's upper x-bound + int iyinl; ///< voxel index of the field's lower y-bound + int iyinu; ///< voxel index of the field's upper y-bound }; -/* Clamp the collimator rectangle to the phantom and find the voxel indices it - covers. A rectangle of zero width in a direction is a pencil beam there. */ +/*! Clamp the collimator rectangle to the phantom and find the voxel indices it + covers. A rectangle of zero width in a direction is a pencil beam there. + + @param source The x/y bounds and ssd must already be filled in; the derived + fields are written by this call. */ void omcSsdSourceInit(struct OmcSsdSource *source); +/*! Run parameters for one cube calculation. */ struct OmcCubeOptions { - int nhist; // total histories - int nbatch; // statistical batches to split them into - int charge; // 0 : photons, -1 : electrons, +1 : positrons + int nhist; ///< total histories + int nbatch; ///< statistical batches to split them into + int charge; ///< 0 : photons, -1 : electrons, +1 : positrons - /* 1 : dose in Gy per incident fluence, 0 : mean deposited energy */ - int outputDose; + int outputDose; ///< 1 : dose in Gy per incident fluence, 0 : mean deposited energy }; +/*! Callbacks omcCalcCube() reports progress through. */ struct OmcCubeCallbacks { - /* About to start a batch; ibatch counts from 0 and firstHistory is the - index of its first history. Optional. Called on the master thread. + /*! About to start a batch. Optional. Called on the master thread. - Return 0 to abandon the calculation: it stops before that batch, tears + @param ibatch Batch index, counting from 0. + @param nbatch Total number of batches. + @param firstHistory Global history index of the batch's first history. + @param user The pointer from struct OmcCubeCallbacks::user, untouched. + @return 0 to abandon the calculation: it stops before that batch, tears its state down and returns 0 without touching dose[] or uncertainty[]. There is no partial result to keep -- the batches are averaged, so a run that stopped halfway would be a dose with no meaning. Return nonzero to carry on. */ int (*batch)(int ibatch, int nbatch, uint64_t firstHistory, void *user); - void *user; + void *user; ///< passed back to the callback, untouched }; -/* What the run did, for hosts that want to report it. Optional. */ +/*! What the run did, for hosts that want to report it. Optional. */ struct OmcCubeSummary { - int nhist; // histories actually run, rounded to whole batches - int nperbatch; - double energyFraction; // deposited energy over incident kinetic energy + int nhist; ///< histories actually run, rounded to whole batches + int nperbatch; ///< histories per batch + double energyFraction; ///< deposited energy over incident kinetic energy }; -/* Transport the histories and write the results into dose[] and, unless it is - NULL, uncertainty[]. Both are supplied by the caller and hold one entry per - voxel, indexed like the phantom: ix + iy*isize + iz*isize*jsize. - - uncertainty is the RELATIVE uncertainty of the dose in that voxel, and is - 0.9999999 wherever nothing was deposited -- the convention the .3ddose format +/*! Transport the histories and write the results into dose[] and, unless it + is NULL, uncertainty[]. Both are supplied by the caller and hold one entry + per voxel, indexed like the phantom: `ix + iy*isize + iz*isize*jsize`. + + @param options Run parameters. + @param source The collimated point source, already passed through + omcSsdSourceInit(). + @param spectrum Source energy spectrum. + @param dose Caller-supplied array of `isize*jsize*ksize` entries. + @param uncertainty Caller-supplied array of the same size, or `NULL`. Holds + the RELATIVE uncertainty of the dose in that voxel, and is 0.9999999 + wherever nothing was deposited -- the convention the .3ddose format expects. - - Returns nonzero when the run finished, 0 when the batch callback stopped it. */ + @param callbacks Progress reporting; see struct OmcCubeCallbacks. + @param summary Optional; filled in with what the run did. + @return Nonzero when the run finished, 0 when the batch callback stopped + it. */ int omcCalcCube(const struct OmcCubeOptions *options, const struct OmcSsdSource *source, const struct OmcSpectrum *spectrum, diff --git a/src/omc_engine_dij.h b/src/omc_engine_dij.h index 1447109..9470521 100644 --- a/src/omc_engine_dij.h +++ b/src/omc_engine_dij.h @@ -19,7 +19,8 @@ along with this program. If not, see . *****************************************************************************/ -/****************************************************************************** +/*! + @file omc_engine_dij - Dose influence matrix for a set of beamlets. This is what the matRad interface calculates: each beamlet is an aperture @@ -40,8 +41,8 @@ the scoring arrays, the random number generators and the particle stacks, which it sets up and tears down per call. - Like the rest of ompMC this is a singleton: one calculation at a time per - process, since the transport state lives in globals. + @warning Like the rest of ompMC this is a singleton: one calculation at a + time per process, since the transport state lives in globals. *****************************************************************************/ #ifndef OMC_ENGINE_DIJ_H @@ -55,47 +56,66 @@ struct OmcSpectrum; +/*! Run parameters for one Dij calculation. */ struct OmcDijOptions { - int nhist; // total histories per beamlet - int nbatch; // statistical batches to split them into - int charge; // 0 : photons, -1 : electrons, +1 : positrons + int nhist; ///< total histories per beamlet + int nbatch; ///< statistical batches to split them into + int charge; ///< 0 : photons, -1 : electrons, +1 : positrons - /* Voxels below this fraction of the beamlet's maximum dose are dropped + /*! Voxels below this fraction of the beamlet's maximum dose are dropped from the column rather than reported. */ double relDoseThreshold; - enum OmcSourceGeometry sourceGeometry; - double sourceGaussianWidth; // standard deviation in cm, GAUSSIAN only + enum OmcSourceGeometry sourceGeometry; ///< POINT or GAUSSIAN, see omc_source_beamlet.h + double sourceGaussianWidth; ///< standard deviation in cm, GAUSSIAN only - int wantVariance; // also report the variance of the mean + int wantVariance; ///< also report the variance of the mean }; +/*! Callbacks omcCalcDij() reports results and progress through. */ struct OmcDijCallbacks { - /* One finished beamlet. voxels holds nvoxels grid indices in ascending - order, 0 based, and dose the dose in Gy in each of them; variance holds - the variance of the mean where it was asked for and is NULL otherwise. - All three arrays belong to the engine and are only valid for the duration - of the call. - - Called on the master thread, outside any parallel region. */ + /*! One finished beamlet. + + @param ibeamlet Index of the beamlet just finished. + @param nvoxels Number of entries in @p voxels, @p dose and @p variance. + @param voxels @p nvoxels grid indices in ascending order, 0 based. + @param dose Dose in Gy, one entry per voxel in @p voxels. + @param variance Variance of the mean, one entry per voxel in @p voxels, + or `NULL` if it was not asked for. + @param user The pointer from struct OmcDijCallbacks::user, untouched. + + All three arrays belong to the engine and are only valid for the + duration of the call. Called on the master thread, outside any parallel + region. */ void (*beamlet)(int ibeamlet, int nvoxels, const int *voxels, const double *dose, const double *variance, void *user); - /* Fraction of the whole calculation finished, in [0,1]. Optional; called - once per batch and once per beamlet, also on the master thread. + /*! Progress report, called once per batch and once per beamlet, also on + the master thread. - Return 0 to abandon the calculation. It stops after the current batch, + @param fraction Fraction of the whole calculation finished, in [0,1]. + @param user The pointer from struct OmcDijCallbacks::user, untouched. + @return 0 to abandon the calculation. It stops after the current batch, tears its state down and returns normally, having reported fewer beamlets than were asked for -- omcCalcDij() tells the caller how many through its return value, and anything already handed to beamlet() stays - valid. Return nonzero to carry on. */ + valid. Return nonzero to carry on. + + Optional: pass `NULL` to skip progress reporting. */ int (*progress)(double fraction, void *user); - void *user; // passed back to both, untouched + void *user; ///< passed back to both callbacks, untouched }; -/* Returns the number of beamlets reported through the beamlet callback, which - is source->nbeamlets unless the progress callback asked to stop early. */ +/*! Run a Dij calculation. + + @param options Run parameters. + @param source The beamlets to calculate a column for. + @param spectrum Source energy spectrum. + @param callbacks Where the results and progress go; see struct + OmcDijCallbacks. + @return The number of beamlets reported through the beamlet callback, which + is `source->nbeamlets` unless the progress callback asked to stop early. */ int omcCalcDij(const struct OmcDijOptions *options, const struct OmcBeamletSource *source, const struct OmcSpectrum *spectrum, diff --git a/src/omc_engine_forward.h b/src/omc_engine_forward.h index ccdc78e..a9774cd 100644 --- a/src/omc_engine_forward.h +++ b/src/omc_engine_forward.h @@ -19,7 +19,8 @@ along with this program. If not, see . *****************************************************************************/ -/****************************************************************************** +/*! + @file omc_engine_forward - Dose from a whole weighted set of beamlets, in one cube. Forward dose for a fluence map. The beamlets and their geometry are the same @@ -36,10 +37,10 @@ so the cost no longer grows with the number of beamlets, and a beamlet that is closed costs nothing at all. - The weights modulate FLUENCE, not spectrum: a weight of 0.02 for a leaf that - transmits 2% starts 2% of the particles it would otherwise have started, with - the unhardened spectrum. Attenuation through the collimator, its scatter and - the beam hardening that comes with it are not modelled. + @warning The weights modulate FLUENCE, not spectrum: a weight of 0.02 for a + leaf that transmits 2% starts 2% of the particles it would otherwise have + started, with the unhardened spectrum. Attenuation through the collimator, + its scatter and the beam hardening that comes with it are not modelled. Before calling omcCalcForward() the host must have @@ -47,8 +48,10 @@ 2. called initMediaData() and initVrt() (ompmc.h) 3. built the source spectrum (omc_spectrum.h) - and afterwards it owns the cleanup of those. Like the rest of ompMC this is a - singleton: one calculation at a time per process. + and afterwards it owns the cleanup of those. + + @warning Like the rest of ompMC this is a singleton: one calculation at a + time per process. *****************************************************************************/ #ifndef OMC_ENGINE_FORWARD_H @@ -58,63 +61,72 @@ struct OmcSpectrum; +/*! Run parameters for one forward calculation. */ struct OmcForwardOptions { - int nhist; // total histories, over all beamlets together - int nbatch; // statistical batches to split them into - int charge; // 0 : photons, -1 : electrons, +1 : positrons + int nhist; ///< total histories, over all beamlets together + int nbatch; ///< statistical batches to split them into + int charge; ///< 0 : photons, -1 : electrons, +1 : positrons - enum OmcSourceGeometry sourceGeometry; - double sourceGaussianWidth; // standard deviation in cm, GAUSSIAN only + enum OmcSourceGeometry sourceGeometry; ///< POINT or GAUSSIAN, see omc_source_beamlet.h + double sourceGaussianWidth; ///< standard deviation in cm, GAUSSIAN only - /* 1 : dose in Gy for the weights given, 0 : mean deposited energy */ - int outputDose; + int outputDose; ///< 1 : dose in Gy for the weights given, 0 : mean deposited energy }; +/*! Callbacks omcCalcForward() reports progress through. */ struct OmcForwardCallbacks { - /* Fraction of the calculation finished, in [0,1]. Optional; called once - per batch on the master thread, outside any parallel region. + /*! Progress report, called once per batch on the master thread, outside + any parallel region. + + @param fraction Fraction of the calculation finished, in [0,1]. + @param user The pointer from struct OmcForwardCallbacks::user, untouched. + @return 0 to abandon the calculation. It stops after the current batch + and returns 0 without touching dose[] or uncertainty[] -- the batches + are averaged, so a run that stopped halfway would be a dose with no + meaning. Return nonzero to carry on. - Return 0 to abandon the calculation. It stops after the current batch and - returns 0 without touching dose[] or uncertainty[] -- the batches are - averaged, so a run that stopped halfway would be a dose with no meaning. - Return nonzero to carry on. */ + Optional: pass `NULL` to skip progress reporting. */ int (*progress)(double fraction, void *user); - void *user; + void *user; ///< passed back to the callback, untouched }; -/* What the run did, for hosts that want to report it. Optional. */ +/*! What the run did, for hosts that want to report it. Optional. */ struct OmcForwardSummary { - int nhist; // histories actually run, rounded to whole batches - int nperbatch; + int nhist; ///< histories actually run, rounded to whole batches + int nperbatch; ///< histories per batch - int nweighted; // beamlets with a weight above zero - int nsampled; // of those, the ones that got any histories + int nweighted; ///< beamlets with a weight above zero + int nsampled; ///< of those, the ones that got any histories - double totalWeight; // sum of the weights asked for - double sampledWeight; // sum over the beamlets that got histories + double totalWeight; ///< sum of the weights asked for + double sampledWeight; ///< sum over the beamlets that got histories - double energyFraction; // deposited energy over incident kinetic energy + double energyFraction; ///< deposited energy over incident kinetic energy }; -/* Transport the histories and write the results into dose[] and, unless it is - NULL, uncertainty[]. Both are supplied by the caller and hold one entry per - voxel, indexed like the phantom: ix + iy*isize + iz*isize*jsize. - - weights holds one finite, non-negative value per beamlet, and their sum must - also be finite. Its scale carries through to the result: doubling every weight - doubles the dose, and the dose returned is the dose for exactly these weights, - so that it can be held against dij*weights. Beamlets are given histories in - proportion to their weight, and one whose share rounds to zero histories - contributes nothing -- the summary - reports how much weight that was, and a warning goes to the host when it is - more than a thousandth of the total. - - uncertainty is the RELATIVE uncertainty of the dose in that voxel, and is - 0.9999999 wherever nothing was deposited, the same convention omcCalcCube() - follows. - - Returns nonzero when the run finished, 0 when the progress callback stopped +/*! Transport the histories and write the results into dose[] and, unless it + is NULL, uncertainty[]. Both are supplied by the caller and hold one entry + per voxel, indexed like the phantom: `ix + iy*isize + iz*isize*jsize`. + + @param options Run parameters. + @param source The beamlets, with the geometry of struct OmcBeamletSource. + @param weights One finite, non-negative value per beamlet, and their sum + must also be finite. Its scale carries through to the result: doubling + every weight doubles the dose, and the dose returned is the dose for + exactly these weights, so that it can be held against `dij*weights`. + Beamlets are given histories in proportion to their weight, and one whose + share rounds to zero histories contributes nothing -- the summary reports + how much weight that was, and a warning goes to the host when it is more + than a thousandth of the total. + @param spectrum Source energy spectrum. + @param dose Caller-supplied array of `isize*jsize*ksize` entries. + @param uncertainty Caller-supplied array of the same size, or `NULL`. Holds + the RELATIVE uncertainty of the dose in that voxel, and is 0.9999999 + wherever nothing was deposited, the same convention omcCalcCube() follows. + @param callbacks Progress reporting; see struct OmcForwardCallbacks. + @param summary Optional; filled in with what the run did. + @return Nonzero when the run finished, 0 when the progress callback stopped it. */ int omcCalcForward(const struct OmcForwardOptions *options, const struct OmcBeamletSource *source, diff --git a/src/omc_geom.h b/src/omc_geom.h index f5b7b90..c68f824 100644 --- a/src/omc_geom.h +++ b/src/omc_geom.h @@ -19,7 +19,8 @@ along with this program. If not, see . *****************************************************************************/ -/****************************************************************************** +/*! + @file omc_geom - The rectilinear voxel phantom every ompMC user code transports in. howfar(), hownear() and regionIndex() are the geometry side of the contract @@ -30,43 +31,53 @@ all see one phantom. What a user code still owns is FILLING the geometry: omc_dosxyz reads an - .egsphant file, the matRad interface takes cubes from MATLAB, and a Python - host will take numpy arrays. Whoever fills it must set every field of struct + .egsphant file, the matRad interface takes cubes from MATLAB, and the Python + host takes numpy arrays. Whoever fills it must set every field of struct Geom below, including the reciprocal spacings, before calling initRegions(). *****************************************************************************/ #ifndef OMC_GEOM_H #define OMC_GEOM_H +/*! The rectilinear voxel phantom. A host fills every field, including the + reciprocal spacings (see omcGeomDetectSpacing()), before calling + initRegions(). */ struct Geom { - int *med_indices; // index of the media in each voxel - double *med_densities; // density of the medium in each voxel - - int isize; // number of voxels on each direction - int jsize; - int ksize; - - double *xbounds; // boundaries of voxels on each direction - double *ybounds; - double *zbounds; - - double dxi, dyi, dzi; /* reciprocal grid spacing per axis when that - axis is uniform, 0.0 when it is not; lets - regionIndex() locate a point with one - multiplication instead of a binary search */ + int *med_indices; ///< index of the medium in each voxel + double *med_densities; ///< density of the medium in each voxel + + int isize; ///< number of voxels along x + int jsize; ///< number of voxels along y + int ksize; ///< number of voxels along z + + double *xbounds; ///< boundaries of voxels along x, isize+1 values + double *ybounds; ///< boundaries of voxels along y, jsize+1 values + double *zbounds; ///< boundaries of voxels along z, ksize+1 values + + /*! Reciprocal grid spacing along x when that axis is uniform, 0.0 + otherwise; lets regionIndex() locate a point with one multiplication + instead of a binary search. Filled by omcGeomDetectSpacing(). */ + double dxi; + double dyi; ///< reciprocal grid spacing along y, see #dxi + double dzi; ///< reciprocal grid spacing along z, see #dxi }; +/*! The phantom every user code fills in and passes to initRegions(). */ extern struct Geom geometry; -/* Fill dxi/dyi/dzi from the bounds already stored in the struct. Every loader - has to do this, and doing it in one place keeps a new one from forgetting and - quietly losing the fast point location. */ +/*! Fill dxi/dyi/dzi from the bounds already stored in the struct. Every + loader has to do this, and doing it in one place keeps a new one from + forgetting and quietly losing the fast point location. */ void omcGeomDetectSpacing(void); -/* Set up the per region transport parameters from the filled geometry: medium - index and density scaling per voxel, per medium cut-offs clamped to what the - PEGS data supports, and the per medium maximum density ratio the Woodcock - majorant needs. Reads the "global ecut" and "global pcut" input items. */ +/*! Set up the per region transport parameters from the filled geometry: + medium index and density scaling per voxel, per medium cut-offs clamped to + what the PEGS data supports, and the per medium maximum density ratio the + Woodcock majorant needs. Reads the "global ecut" and "global pcut" input + items. + + @pre The global #geometry is completely filled in, including a call to + omcGeomDetectSpacing(). */ void initRegions(void); #endif diff --git a/src/omc_host.h b/src/omc_host.h index ea3c8a6..0533575 100644 --- a/src/omc_host.h +++ b/src/omc_host.h @@ -19,7 +19,8 @@ along with this program. If not, see . *****************************************************************************/ -/****************************************************************************** +/*! + @file omc_host - How shared ompMC code talks back to whoever is embedding it. Code that is meant to be used from more than one host -- a command line @@ -27,18 +28,19 @@ mexErrMsgIdAndTxt() directly. It calls omcLog() and omcFail() instead, and the host installs the sinks that give those meaning. - THREADING: both sinks are called on the master thread only, never from inside - an OpenMP parallel region. Hosts may therefore call back into a runtime that - has no business being entered from a worker thread -- MATLAB's mexPrintf, or a - Python callable under the GIL. Diagnostics that transport code emits from - inside a parallel region keep going straight to stdout instead, which is why - ompmc.c still uses printf() there. + @warning THREADING: both sinks are called on the master thread only, never + from inside an OpenMP parallel region. Hosts may therefore call back into a + runtime that has no business being entered from a worker thread -- MATLAB's + mexPrintf, or a Python callable under the GIL. Diagnostics that transport + code emits from inside a parallel region keep going straight to stdout + instead, which is why ompmc.c still uses printf() there. *****************************************************************************/ #ifndef OMC_HOST_H #define OMC_HOST_H -/* omcFail() never comes back, and saying so lets callers end a function with +/*! @cond OMC_INTERNAL + omcFail() never comes back, and saying so lets callers end a function with it the way they used to end one with exit(), without the compiler asking for a return value it will never need. The format attribute keeps the printf style arguments checked, which the direct printf() calls got for free. */ @@ -53,45 +55,64 @@ #define OMC_NORETURN #define OMC_PRINTF_LIKE(fmtArg, firstArg) #endif +/*! @endcond */ -/* Severity of a message passed to omcLog(). The sink decides what to do with - each level; nothing is filtered on the way there, so that a host can be as - chatty or as quiet as it likes without the shared code knowing. */ +/*! Severity of a message passed to omcLog(). The sink decides what to do + with each level; nothing is filtered on the way there, so that a host can be + as chatty or as quiet as it likes without the shared code knowing. */ enum OmcLogLevel { OMC_LOG_WARNING = 0, - OMC_LOG_INFO, /* progress and summaries */ - OMC_LOG_DETAIL, /* details a curious user might want */ - OMC_LOG_DEBUG /* dumps only useful when something is wrong */ + OMC_LOG_INFO, /**< progress and summaries */ + OMC_LOG_DETAIL, /**< details a curious user might want */ + OMC_LOG_DEBUG /**< dumps only useful when something is wrong */ }; +/*! The sinks a host installs with omcSetHost(). */ struct OmcHost { - /* Receives an already formatted message, without a trailing newline. */ + /*! Receives an already formatted message, without a trailing newline. + + @param level One of enum OmcLogLevel. + @param message The formatted message. + @param user The pointer from struct OmcHost::user, untouched. */ void (*log)(int level, const char *message, void *user); - /* Reports a fatal condition. MUST NOT RETURN: the shared code calls this - where it has no way to carry on, and simply continues into undefined - state if the call comes back. Hosts end it by exiting the process - (command line), throwing out of the call (MATLAB, Octave) or jumping back - to the entry point with longjmp() (Python). omcFail() calls abort() if a - sink returns anyway, so the mistake is loud rather than silent. + /*! Reports a fatal condition. + + @param id Dotted identifier for hosts that can carry one, e.g. + `"ompMC:geometry:badMaterialIndex"`. Hosts that cannot may ignore it. + @param message The formatted message. + @param user The pointer from struct OmcHost::user, untouched. - id is a dotted identifier for hosts that can carry one, e.g. - "ompMC:geometry:badMaterialIndex". Hosts that cannot may ignore it. */ + @warning MUST NOT RETURN: the shared code calls this where it has no way + to carry on, and simply continues into undefined state if the call comes + back. Hosts end it by exiting the process (command line), throwing out + of the call (MATLAB, Octave) or jumping back to the entry point with + longjmp() (Python). omcFail() calls abort() if a sink returns anyway, so + the mistake is loud rather than silent. */ void (*fail)(const char *id, const char *message, void *user); - void *user; /* passed back to both sinks untouched */ + void *user; ///< passed back to both sinks untouched }; -/* Install the sinks. Passing NULL restores the built-in default, which prints - to stdout/stderr and exits the process on a failure -- what a plain command +/*! Install the sinks. + + @param host The sinks to install; the pointer is not retained, the struct + is copied. Passing `NULL` restores the built-in default, which prints to + stdout/stderr and exits the process on a failure -- what a plain command line program wants, and a safe fallback for a host that forgets to install - its own. The pointer is not retained; the struct is copied. */ + its own. */ void omcSetHost(const struct OmcHost *host); -/* Format and hand a message to the log sink. */ +/*! Format and hand a message to the log sink. + + @param level One of enum OmcLogLevel. + @param fmt printf style format string, followed by its arguments. */ void omcLog(int level, const char *fmt, ...) OMC_PRINTF_LIKE(2, 3); -/* Format and hand a fatal message to the fail sink. Does not return. */ +/*! Format and hand a fatal message to the fail sink. Does not return. + + @param id Dotted identifier, see struct OmcHost::fail. + @param fmt printf style format string, followed by its arguments. */ OMC_NORETURN void omcFail(const char *id, const char *fmt, ...) OMC_PRINTF_LIKE(2, 3); diff --git a/src/omc_random.h b/src/omc_random.h index bc431fe..73b6c32 100644 --- a/src/omc_random.h +++ b/src/omc_random.h @@ -21,39 +21,43 @@ along with this program. If not, see . *****************************************************************************/ -/******************************************************************************* -* Counter-based random number generator built on Philox4x32-10 (Salmon, -* Moraes, Dror and Shaw, "Parallel random numbers: as easy as 1, 2, 3", -* SC'11). It replaces the RANMAR port used previously. -* -* The generator is a pure function of a 64 bit key and a 128 bit counter. -* The key comes from the 'rng seeds' input. The high 64 bits of the counter -* hold the global history index, set through setRandomHistory() at the start -* of every particle history; the low 64 bits count the draws within the -* history. Every history therefore owns its own stream of 2^64 numbers, -* determined only by the seeds and the history index -- never by the thread -* that happens to simulate it or by how histories are scheduled. -* -* Before using the RNG, it is needed to initialize the RNG by a call to -* initRandom(). -*******************************************************************************/ +/*! + @file + Counter-based random number generator built on Philox4x32-10 (Salmon, + Moraes, Dror and Shaw, "Parallel random numbers: as easy as 1, 2, 3", + SC'11). It replaces the RANMAR port used previously. + + The generator is a pure function of a 64 bit key and a 128 bit counter. + The key comes from the 'rng seeds' input. The high 64 bits of the counter + hold the global history index, set through setRandomHistory() at the start + of every particle history; the low 64 bits count the draws within the + history. Every history therefore owns its own stream of 2^64 numbers, + determined only by the seeds and the history index -- never by the thread + that happens to simulate it or by how histories are scheduled. + + @warning Before using the RNG, it is needed to initialize the RNG by a call + to initRandom(). +*****************************************************************************/ #include -#define BUFF_SIZE 256 +#define BUFF_SIZE 256 ///< size of the scratch buffer initRandom() reads seeds into -/* Scale factor turning 32 bit words into reals. Exact in binary floating +/*! Scale factor turning 32 bit words into reals. Exact in binary floating point. */ #define TWOM32 (1.0/4294967296.0) +/*! Per-thread generator state: the Philox4x32-10 key and counter, and a + small buffer of already-converted reals. */ struct Random { - uint32_t key[2]; /* base key, taken from the 'rng seeds' input */ - uint32_t ctr[4]; /* ctr[2],ctr[3] hold the history index; ctr[0], + uint32_t key[2]; /**< base key, taken from the 'rng seeds' input */ + uint32_t ctr[4]; /**< ctr[2],ctr[3] hold the history index; ctr[0], ctr[1] count the blocks drawn within the history */ - int buf_pos; /* next unread entry of buf; 4 means empty */ - double buf[4]; /* one Philox block converted to reals in (0,1) */ + int buf_pos; /**< next unread entry of buf; 4 means empty */ + double buf[4]; /**< one Philox block converted to reals in (0,1) */ }; +/*! Per-thread generator state. */ #if defined(_MSC_VER) /* use __declspec(thread) instead of threadprivate to avoid error C3053. More information in: @@ -64,31 +68,45 @@ struct Random { #pragma omp threadprivate(rng) #endif #ifndef M_PI - #define M_PI 3.14159265358979323846 + #define M_PI 3.14159265358979323846 ///< pi, for compilers whose math.h omits it #endif -/* Read the 'rng seeds' input into the thread-local key and leave the +/*! Read the 'rng seeds' input into the thread-local key and leave the generator on a sentinel stream no real history uses. Call once per thread before any setRandom(). */ void initRandom(void); -/* Point the generator at the stream owned by global history index ihist. +/*! Point the generator at the stream owned by global history index ihist. Call at the start of every particle history; the index must be unique over - the whole run (across batches, and beamlets where applicable). */ + the whole run (across batches, and beamlets where applicable). + + @param ihist Global history index. */ void setRandomHistory(uint64_t ihist); -/* Get a single floating random number in (0,1) from the current stream */ +/*! @return A single floating random number in (0,1) from the current + stream. */ double setRandom(void); -/* One Philox4x32-10 block: 128 bit counter and 64 bit key in, four 32 bit - words out. Exposed for verification against the published test vectors. */ +/*! One Philox4x32-10 block: 128 bit counter and 64 bit key in, four 32 bit + words out. Exposed for verification against the published test vectors. + + @param ctr 128 bit counter, as four 32 bit words. + @param key 64 bit key, as two 32 bit words. + @param out Four 32 bit output words. */ void philox4x32(const uint32_t ctr[4], const uint32_t key[2], uint32_t out[4]); +/*! Release any resources initRandom() allocated. */ void cleanRandom(void); +/*! @return A normally distributed random number with mean mu and standard + deviation sigma. */ double setStandardNormalRandom(const double mu, const double sigma); +/*! Box-Muller transform: two independent standard normal deviates from two + uniform ones. + + @param rndnormal Filled with the two deviates. */ void boxMuller(double rndnormal[2]); /******************************************************************************/ diff --git a/src/omc_score.h b/src/omc_score.h index a5fd6ad..22515ad 100644 --- a/src/omc_score.h +++ b/src/omc_score.h @@ -21,44 +21,47 @@ along with this program. If not, see . *****************************************************************************/ -/******************************************************************************* -* Energy scoring, shared by the user codes. -* -* The dose from one beamlet occupies a small fraction of the grid, often well -* under one percent, so anything that sweeps the whole grid once per batch or -* once per beamlet costs far more than the scoring itself -- and a Dij run -* does exactly that, nbeamlets times over. Every voxel that receives energy is -* therefore recorded in a list, and the accumulation, output and reset steps -* walk that list instead of the grid. -* -* ausgab() itself keeps the plain atomic add. Combining deposits per voxel in -* thread local state first was tried and measured slower: an electron does not -* stay in one voxel for long enough runs of steps to pay for touching thread -* local state on every call, and the extra branch and TLS access cost about -* ten percent on omc_dosxyz while saving nothing measurable. -* -* Call order over a beamlet: -* -* initScore(gridsize) once -* ... per batch: -* parallel { initHistory(); shower(); } -* accumEndep(scale) -* scoreBeamVoxels(&list) to emit the beamlet's results -* resetBeamScore() before the next beamlet -* cleanScore() once -*******************************************************************************/ +/*! + @file + Energy scoring, shared by the user codes. + + The dose from one beamlet occupies a small fraction of the grid, often well + under one percent, so anything that sweeps the whole grid once per batch or + once per beamlet costs far more than the scoring itself -- and a Dij run + does exactly that, nbeamlets times over. Every voxel that receives energy is + therefore recorded in a list, and the accumulation, output and reset steps + walk that list instead of the grid. + + ausgab() itself keeps the plain atomic add. Combining deposits per voxel in + thread local state first was tried and measured slower: an electron does not + stay in one voxel for long enough runs of steps to pay for touching thread + local state on every call, and the extra branch and TLS access cost about + ten percent on omc_dosxyz while saving nothing measurable. + + Call order over a beamlet: + + initScore(gridsize) once + ... per batch: + parallel { initHistory(); shower(); } + accumEndep(scale) + scoreBeamVoxels(&list) to emit the beamlet's results + resetBeamScore() before the next beamlet + cleanScore() once +*****************************************************************************/ +/*! Accumulated dose and the voxel list it was accumulated over. */ struct Score { - double ensrc; // total energy from source + double ensrc; ///< total energy from source - int gridsize; // voxel count; regions run 1..gridsize, with - // region 0 reserved for outside the geometry + /*! voxel count; regions run 1..gridsize, with region 0 reserved for + outside the geometry */ + int gridsize; - double *endep; // energy deposited during the current batch - double *accum_endep; // accumulated across batches - double *accum_endep2; // accumulated squares, for the batch variance + double *endep; ///< energy deposited during the current batch + double *accum_endep; ///< accumulated across batches + double *accum_endep2; ///< accumulated squares, for the batch variance - /* Voxels this beamlet has deposited energy in, i.e. since the last + /*! Voxels this beamlet has deposited energy in, i.e. since the last resetBeamScore(). The flag keeps the list duplicate free; see omc_score.c for how it stays correct under concurrent updates. @@ -68,43 +71,57 @@ struct Score { save nothing in the accumulation while making every batch re-enter the lock once for every voxel it hits. */ int *beam_list; - int beam_count; - int beam_sorted; // beam_list is already in ascending order - unsigned char *beam_flag; + int beam_count; ///< number of entries in use in beam_list + int beam_sorted; ///< beam_list is already in ascending order + unsigned char *beam_flag; ///< per-voxel "already in beam_list" flag }; +/*! The (single) accumulated score for the calculation in progress. */ extern struct Score score; -/* gridsize is the number of voxels, not counting region 0. */ +/*! @param gridsize The number of voxels, not counting region 0. */ void initScore(int gridsize); + +/*! Release the accumulators initScore() allocated. */ void cleanScore(void); -/* Add to the source energy tally. Safe to call from inside a parallel - region, unlike a bare score.ensrc += ein. */ +/*! Add to the source energy tally. Safe to call from inside a parallel + region, unlike a bare `score.ensrc += ein`. + + @param ein Energy to add, in MeV. */ void scoreSource(double ein); -/* Fold the current batch into the accumulators, multiplying by scale, then +/*! Fold the current batch into the accumulators, multiplying by scale, then clear it. Costs O(voxels this beamlet has touched), not O(gridsize). Call - outside any parallel region. */ + outside any parallel region. + + @param scale Factor the current batch's deposits are multiplied by before + being added to the accumulators. */ void accumEndep(double scale); -/* The voxels this beamlet has deposited energy in, ascending so that callers - can emit them straight into a column of a CSC sparse matrix. Returns the - count and, through list, the indices. */ +/*! The voxels this beamlet has deposited energy in, ascending so that + callers can emit them straight into a column of a CSC sparse matrix. + + @param list Set to point at the (engine-owned) list of voxel indices. + @return The number of entries in @p list. */ int scoreBeamVoxels(const int **list); -/* Zero the accumulators over this beamlet's touched set and begin a new +/*! Zero the accumulators over this beamlet's touched set and begin a new beamlet. Costs O(voxels touched). */ void resetBeamScore(void); -/* Turn what the batches accumulated into a dense dose cube and its relative +/*! Turn what the batches accumulated into a dense dose cube and its relative uncertainty, one entry per voxel, indexed like the phantom: - ix + iy*isize + iz*isize*jsize. uncertainty may be NULL. - - outputDose selects Gy (1) or mean deposited energy (0), and incFluence is - what the accumulated energy is divided by -- see the comment on the - definition. Walks the whole grid, so that empty voxels come out with the - 0.9999999 the .3ddose format expects. Call outside any parallel region. */ + `ix + iy*isize + iz*isize*jsize`. Walks the whole grid, so that empty + voxels come out with the 0.9999999 the .3ddose format expects. Call outside + any parallel region. + + @param nbatch Number of batches accumulated. + @param incFluence What the accumulated energy is divided by -- see the + comment on the definition. + @param outputDose Selects Gy (1) or mean deposited energy (0). + @param dose Caller-supplied array of gridsize entries. + @param uncertainty Caller-supplied array of gridsize entries, or `NULL`. */ void omcScoreToCube(int nbatch, double incFluence, int outputDose, double *dose, double *uncertainty); diff --git a/src/omc_source_beamlet.h b/src/omc_source_beamlet.h index 3919739..7b9d89e 100644 --- a/src/omc_source_beamlet.h +++ b/src/omc_source_beamlet.h @@ -19,7 +19,8 @@ along with this program. If not, see . *****************************************************************************/ -/****************************************************************************** +/*! + @file omc_source_beamlet - Particles starting on a beamlet aperture. The source model the matRad interface uses: a beam has a source point, a @@ -42,59 +43,64 @@ struct OmcSpectrum; -/* Where on the source the particles start. POINT is the classic point source; - GAUSSIAN spreads the starting point over the collimator plane, which softens - the penumbra. */ +/*! Where on the source the particles start. */ enum OmcSourceGeometry { - OMC_SOURCE_POINT = 0, - OMC_SOURCE_GAUSSIAN + OMC_SOURCE_POINT = 0, /**< the classic point source */ + OMC_SOURCE_GAUSSIAN /**< spreads the starting point over the + collimator plane, which softens the penumbra */ }; -/* The beamlets. Per beam: the source position. Per beamlet: which beam it +/*! The beamlets. Per beam: the source position. Per beamlet: which beam it belongs to, and the corner plus two edge vectors of its aperture rectangle at - isocentre. All arrays belong to the caller and must outlive the call. */ + isocentre. + + @warning All arrays belong to the caller and must outlive the call. */ struct OmcBeamletSource { - int nbeamlets; - const int *ibeam; // index of the beam of each beamlet, 0 based + int nbeamlets; ///< number of beamlets + const int *ibeam; ///< index of the beam of each beamlet, 0 based - const double *xsource; // coordinates of the source of each beam - const double *ysource; - const double *zsource; + const double *xsource; ///< x coordinate of the source of each beam + const double *ysource; ///< y coordinate of the source of each beam + const double *zsource; ///< z coordinate of the source of each beam - const double *xcorner; // coordinates of the beamlet corner - const double *ycorner; - const double *zcorner; + const double *xcorner; ///< x coordinate of the beamlet corner + const double *ycorner; ///< y coordinate of the beamlet corner + const double *zcorner; ///< z coordinate of the beamlet corner - const double *xside1; // first edge vector of the beamlet - const double *yside1; - const double *zside1; + const double *xside1; ///< x component of the first edge vector of the beamlet + const double *yside1; ///< y component of the first edge vector of the beamlet + const double *zside1; ///< z component of the first edge vector of the beamlet - const double *xside2; // second edge vector of the beamlet - const double *yside2; - const double *zside2; + const double *xside2; ///< x component of the second edge vector of the beamlet + const double *yside2; ///< y component of the second edge vector of the beamlet + const double *zside2; ///< z component of the second edge vector of the beamlet }; -/* Everything the sampling needs that does not change from history to history. - An engine fills this once, before its first parallel region, and hands it to - omcBeamletSample() unchanged from then on. */ +/*! Everything the sampling needs that does not change from history to + history. An engine fills this once, before its first parallel region, and + hands it to omcBeamletSample() unchanged from then on. */ struct OmcBeamletSampler { - const struct OmcBeamletSource *source; - const struct OmcSpectrum *spectrum; + const struct OmcBeamletSource *source; ///< the beamlets + const struct OmcSpectrum *spectrum; ///< source energy spectrum - int charge; // 0 : photons, -1 : electrons, +1 : positrons - enum OmcSourceGeometry geometry; - double gaussianWidth; // standard deviation in cm, GAUSSIAN only + int charge; ///< 0 : photons, -1 : electrons, +1 : positrons + enum OmcSourceGeometry geometry; ///< POINT or GAUSSIAN + double gaussianWidth; ///< standard deviation in cm, GAUSSIAN only }; -/* Put one primary particle of beamlet ibeamlet on the (thread local) stack, +/*! Put one primary particle of beamlet ibeamlet on the (thread local) stack, already transported to the phantom surface and with its region index found. - weight becomes the particle's statistical weight, and also scales what the - history contributes to the incident energy tally. Pass 1.0 for an unweighted - history; omc_engine_forward passes the beamlet's share of the fluence. + @param sampler Sampling parameters, unchanged since the caller filled them. + @param ibeamlet Index of the beamlet to sample from. + @param weight Becomes the particle's statistical weight, and also scales + what the history contributes to the incident energy tally. Pass 1.0 for an + unweighted history; omc_engine_forward passes the beamlet's share of the + fluence. - Runs inside the parallel history loop, so it touches nothing but the thread's - own stack, its own random number generator, and the read-only sampler. */ + @warning Runs inside the parallel history loop, so it touches nothing but + the thread's own stack, its own random number generator, and the read-only + sampler. */ void omcBeamletSample(const struct OmcBeamletSampler *sampler, int ibeamlet, double weight); diff --git a/src/omc_spectrum.h b/src/omc_spectrum.h index e1ca44d..0316e7c 100644 --- a/src/omc_spectrum.h +++ b/src/omc_spectrum.h @@ -19,7 +19,8 @@ along with this program. If not, see . *****************************************************************************/ -/****************************************************************************** +/*! + @file omc_spectrum - The energy distribution of the source particles. A spectrum is a histogram: nbins bins, bin i running from the upper energy of @@ -40,46 +41,69 @@ #include -/* What counts[] means. Counts per MeV are converted to counts per bin on the - way in, by scaling with the bin widths. */ +/*! What counts[] means when building a spectrum from a histogram. Counts + per MeV are converted to counts per bin on the way in, by scaling with the + bin widths. */ enum OmcSpectrumMode { OMC_SPECTRUM_COUNTS_PER_BIN = 0, OMC_SPECTRUM_COUNTS_PER_MEV = 1 }; +/*! A source energy spectrum, ready to sample from. */ struct OmcSpectrum { - int monoenergetic; /* 1 : every particle starts at energy */ - double energy; /* the energy, when monoenergetic */ + int monoenergetic; ///< 1 : every particle starts at #energy + double energy; ///< the energy, in MeV, when #monoenergetic - double deltak; /* number of elements in the inverse CDF */ - double *cdfinv1; /* lower energy of the bin an element falls in */ - double *cdfinv2; /* width of that bin */ + double deltak; ///< number of elements in the inverse CDF + double *cdfinv1; ///< lower energy of the bin an element falls in + double *cdfinv2; ///< width of that bin }; /* All of these leave the spectrum ready to sample from, and take ownership of nothing: the arrays passed in may be freed by the caller afterwards. */ +/*! @param spectrum Filled in as a monoenergetic spectrum. + @param energy The energy, in MeV, every particle starts at. */ void omcSpectrumMonoenergetic(struct OmcSpectrum *spectrum, double energy); -/* nbins bins with the given upper energies, ascending and all above emin, and - non-negative counts that sum to something positive. A caller that cannot - guarantee that should check first -- the failures here are reported through - omcFail(), which does not return. */ +/*! Build the inverse-CDF sampling tables from a histogram. + + @param spectrum Filled in from the histogram. + @param upperEnergy Upper energy of each bin, ascending and all above @p + emin. + @param counts Count (or count density, see @p mode) of each bin, + non-negative and summing to something positive. + @param nbins Number of bins, i.e. the length of @p upperEnergy and @p + counts. + @param emin Lower edge of the first bin, in MeV. + @param mode One of enum OmcSpectrumMode. + + @warning A caller that cannot guarantee the constraints above should check + first -- the failures here are reported through omcFail(), which does not + return. */ void omcSpectrumFromHistogram(struct OmcSpectrum *spectrum, const double *upperEnergy, const double *counts, int nbins, double emin, int mode); -/* Read an EGSnrc style .spectrum file: a title line, then "nbins emin mode", - then one "upperEnergy count" pair per line. */ +/*! Read an EGSnrc style .spectrum file: a title line, then "nbins emin mode", + then one "upperEnergy count" pair per line. + + @param spectrum Filled in from the file. + @param path Path to the .spectrum file. */ void omcSpectrumFromFile(struct OmcSpectrum *spectrum, const char *path); +/*! Release the sampling tables omcSpectrumFromHistogram() or + omcSpectrumFromFile() allocated. */ void omcSpectrumFree(struct OmcSpectrum *spectrum); -/* Sample a kinetic energy in MeV. Called once per history, so it is inline +/*! Sample a kinetic energy in MeV. Called once per history, so it is inline rather than a call into another object file; the arithmetic is unchanged from when it sat in the user codes. - It draws its own random numbers, and deliberately draws NONE for a + @param spectrum The spectrum to sample from. + @return A kinetic energy in MeV. + + @warning It draws its own random numbers, and deliberately draws NONE for a monoenergetic source. That is not just an optimization: the random stream is indexed per history, so drawing two numbers that are then thrown away would shift every later draw in the history and change the result of an otherwise diff --git a/src/omc_utilities.h b/src/omc_utilities.h index a7692f7..215c122 100644 --- a/src/omc_utilities.h +++ b/src/omc_utilities.h @@ -3,7 +3,7 @@ /****************************************************************************** ompMC - An OpenMP parallel implementation for Monte Carlo particle transport simulations - + Copyright (C) 2020 Edgardo Doerner (edoerner@fis.puc.cl) @@ -21,88 +21,112 @@ along with this program. If not, see . *****************************************************************************/ -/******************************************************************************/ -/* Timing utilities. If OpenMP is enabled it calculates the wall time through - omp_get_wtime() function. Otherwise, it calculates CPU time through the clock() - function, available in time.h library. */ +/*! + @file + Small utilities shared by the user codes: timing, input-file parsing, and + voxel geometry helpers. +*****************************************************************************/ +/*! @return The wall time in seconds if OpenMP is enabled (via + `omp_get_wtime()`), otherwise the CPU time (via `clock()`, from `time.h`). */ double omc_get_time(); -/******************************************************************************/ -/******************************************************************************/ -/* A simple C/C++ class to parse input files and return requested key value -https://github.com/bmaynard/iniReader */ +/*! + Input-file parsing: a simple C/C++ class to parse input files and return + requested key value pairs. https://github.com/bmaynard/iniReader +*/ -#define BUFFER_SIZE 256 -#define INPUT_PAIRS 80 -#define INPUT_EXT ".inp" // extension of input files +#define BUFFER_SIZE 256 ///< maximum length of one input key or value, including the terminator +#define INPUT_PAIRS 80 ///< maximum number of key/value pairs #input_items holds +#define INPUT_EXT ".inp" ///< extension of input files -/* A path assembled by appending one of the data file names to a folder read - from the input table. The folder is a value, so it is at most BUFFER_SIZE-1 - characters, and the names appended to it are short. Sizing these buffers by - hand is how a Python package installed under a long temporary directory used - to run off the end of a 128 byte array. */ +/*! A path assembled by appending one of the data file names to a folder read + from the input table. The folder is a value, so it is at most + `BUFFER_SIZE-1` characters, and the names appended to it are short. Sizing + these buffers by hand is how a Python package installed under a long + temporary directory used to run off the end of a 128 byte array. */ #define PATH_SIZE (BUFFER_SIZE + 32) -/* Parse a configuration file */ +/*! Parse a configuration file into #input_items. + + @param file_name Path to the input file. */ void parseInputFile(char *file_name); -/* Copy the value of the selected input item to the char pointer */ +/*! Copy the value of the selected input item to the char pointer. + + @param dest Buffer the value is copied into; must be at least + `BUFFER_SIZE` bytes. + @param key The key to look up. + @return Nonzero if @p key was found. */ int getInputValue(char *dest, char *key); -/* Set one key/value pair directly, for hosts that get their configuration - from somewhere other than a file -- a MATLAB struct, a Python dict. Replaces - the value when the key is already there. Both strings are copied, and are - truncated at BUFFER_SIZE-1 characters. */ +/*! Set one key/value pair directly, for hosts that get their configuration + from somewhere other than a file -- a MATLAB struct, a Python dict. + Replaces the value when the key is already there. + + @param key The key, copied and truncated at `BUFFER_SIZE-1` characters. + @param value The value, copied and truncated at `BUFFER_SIZE-1` characters. */ void omcSetInputValue(const char *key, const char *value); -/* Forget every key/value pair. Hosts that stay resident between runs -- a MEX - file, a Python module -- have to start each run from a clean table rather - than inheriting the previous one. */ +/*! Forget every key/value pair. Hosts that stay resident between runs -- a + MEX file, a Python module -- have to start each run from a clean table + rather than inheriting the previous one. */ void omcClearInputValues(void); -/* Returns nonzero if line is a string containing only whitespace or is empty */ +/*! @param line The string to test. + @return Nonzero if @p line contains only whitespace or is empty. */ int lineBlack(char *line); -/* Remove white spaces from string str_untrimmed and saves the results in - str_trimmed. Useful for string input values, such as file names */ +/*! Remove white spaces from a string. Useful for string input values, such + as file names. + + @param str_trimmed Destination buffer for the result. + @param str_untrimmed The string to remove whitespace from. */ void removeSpaces(char* str_trimmed, const char* str_untrimmed); +/*! One key/value pair of #input_items. */ struct inputItems { - char key[BUFFER_SIZE]; - char value[BUFFER_SIZE]; + char key[BUFFER_SIZE]; ///< the key, NUL terminated + char value[BUFFER_SIZE]; ///< its value, NUL terminated }; -/* The key/value table itself. Declared here rather than left for each user +/*! The key/value table itself. Declared here rather than left for each user code to declare extern for itself, because that is how the two halves of the invariant below drifted apart in the first place. - input_idx is the NUMBER of pairs stored, and they occupy input_items[0] up to - input_items[input_idx - 1]. An empty table is input_idx == 0, with no slot to - look at -- which is what makes "is this table empty" answerable at all. + #input_idx is the NUMBER of pairs stored, and they occupy `input_items[0]` + up to `input_items[input_idx - 1]`. An empty table is `input_idx == 0`, + with no slot to look at -- which is what makes "is this table empty" + answerable at all. - Anything filling the table directly rather than through omcSetInputValue() - has to leave it that way. */ + @warning Anything filling the table directly rather than through + omcSetInputValue() has to leave it that way. */ extern struct inputItems input_items[INPUT_PAIRS]; -extern int input_idx; +extern int input_idx; ///< number of pairs stored in #input_items +/*! + Voxel geometry helpers, shared by the user codes. Both omcDecodeRegion() and + omcFindVoxelIndex() are on the transport hot path -- omcDecodeRegion() runs + on every howfar() and hownear() call -- so they are inline in the header + rather than a call into another translation unit. Keeping them here also + makes them reachable from the unit tests. +*/ -/******************************************************************************/ - -/******************************************************************************/ -/* Voxel geometry helpers shared by the user codes. Both are on the transport - hot path -- omcDecodeRegion() runs on every howfar() and hownear() call -- so - they are inline in the header rather than a call into another translation - unit. Keeping them here also makes them reachable from the unit tests. */ - -/* Decode a region number into its voxel indices along each axis. Regions are - numbered 1 + ix + iy*imax + iz*imax*jmax, with 0 reserved for "outside the - geometry"; irl must be >= 1. +/*! Decode a region number into its voxel indices along each axis. Regions are + numbered `1 + ix + iy*imax + iz*imax*jmax`, with 0 reserved for "outside the + geometry". Written with two integer divisions rather than the three the arithmetic - suggests: the quotient of the first division is imax*(iy + iz*jmax)/imax, + suggests: the quotient of the first division is `imax*(iy + iz*jmax)/imax`, i.e. exactly the combined y,z index, so the second division can work on that - directly. Each division pairs with its own remainder into one instruction. */ + directly. Each division pairs with its own remainder into one instruction. + + @param irl Region number, must be >= 1. + @param imax Number of voxels along x. + @param jmax Number of voxels along y. + @param ix Set to the voxel index along x. + @param iy Set to the voxel index along y. + @param iz Set to the voxel index along z. */ static inline void omcDecodeRegion(int irl, int imax, int jmax, int *ix, int *iy, int *iz) { @@ -114,10 +138,14 @@ static inline void omcDecodeRegion(int irl, int imax, int jmax, *iy = irxy - (*iz)*jmax; } -/* Index of the voxel along one axis containing pos, i.e. the smallest i in - [0, n-1] with bounds[i+1] >= pos. bounds holds n+1 ascending values. - Positions outside the grid clamp to the first or last voxel rather than - running off the end of bounds[]. */ +/*! Index of the voxel along one axis containing pos, i.e. the smallest i in + `[0, n-1]` with `bounds[i+1] >= pos`. Positions outside the grid clamp to + the first or last voxel rather than running off the end of @p bounds. + + @param bounds `n+1` ascending values. + @param n Number of voxels along this axis. + @param pos Position along this axis. + @return The voxel index. */ static inline int omcFindVoxelIndex(const double *bounds, int n, double pos) { int lo = 0; @@ -136,10 +164,13 @@ static inline int omcFindVoxelIndex(const double *bounds, int n, double pos) { return lo; } -/* Reciprocal of the grid spacing if the n+1 values in bounds are uniformly - spaced to within a relative tolerance, 0.0 otherwise. Evaluated once at - initialization; the tolerance absorbs the single-precision noise that - phantom files carry in their boundary lists. */ +/*! Evaluated once at initialization; the tolerance absorbs the + single-precision noise that phantom files carry in their boundary lists. + + @param bounds `n+1` ascending values. + @param n Number of voxels along this axis. + @return The reciprocal of the grid spacing if the values in @p bounds are + uniformly spaced to within a relative tolerance, 0.0 otherwise. */ static inline double omcUniformSpacingInv(const double *bounds, int n) { double dx = (bounds[n] - bounds[0])/(double)n; @@ -154,10 +185,17 @@ static inline double omcUniformSpacingInv(const double *bounds, int n) { return 1.0/dx; } -/* Voxel index of pos along one axis: a single multiplication on a uniform - grid (invdx from omcUniformSpacingInv()), the binary search otherwise. The - clamp keeps in-range results for positions on the outer boundaries; callers - reject positions outside the grid before asking. */ +/*! Voxel index of pos along one axis: a single multiplication on a uniform + grid (@p invdx from omcUniformSpacingInv()), the binary search otherwise. + The clamp keeps in-range results for positions on the outer boundaries; + callers reject positions outside the grid before asking. + + @param bounds `n+1` ascending values. + @param n Number of voxels along this axis. + @param invdx Reciprocal grid spacing from omcUniformSpacingInv(), or 0.0 for + a non-uniform grid. + @param pos Position along this axis. + @return The voxel index. */ static inline int omcVoxelIndexFast(const double *bounds, int n, double invdx, double pos) { @@ -175,7 +213,7 @@ static inline int omcVoxelIndexFast(const double *bounds, int n, return omcFindVoxelIndex(bounds, n, pos); } -/* Thread-local memo used by the user codes' howfar()/hownear() to keep the +/*! Thread-local memo used by the user codes' howfar()/hownear() to keep the integer divisions of omcDecodeRegion() and the floating point divisions by the direction cosines off the per-step hot path. @@ -196,13 +234,20 @@ static inline int omcVoxelIndexFast(const double *bounds, int n, and is rejected by the callers before any lookup, and no transported particle has the zero direction. */ struct OmcGeomCache { - int irl; /* region the indices below belong to; 0 = empty */ - int ix, iy, iz; - - double u, v, w; /* direction the reciprocals below belong to */ - double ui, vi, wi; + int irl; /**< region the indices below belong to; 0 = empty */ + int ix; /**< voxel index along x of region #irl */ + int iy; /**< voxel index along y of region #irl */ + int iz; /**< voxel index along z of region #irl */ + + double u; /**< x component of the direction the reciprocals below belong to */ + double v; /**< y component of the direction the reciprocals below belong to */ + double w; /**< z component of the direction the reciprocals below belong to */ + double ui; /**< reciprocal of #u */ + double vi; /**< reciprocal of #v */ + double wi; /**< reciprocal of #w */ }; +/*! Per-thread instance of struct OmcGeomCache. */ #if defined(_MSC_VER) extern __declspec(thread) struct OmcGeomCache omc_geom_cache; #else @@ -210,6 +255,8 @@ struct OmcGeomCache { #pragma omp threadprivate(omc_geom_cache) #endif +/*! omcDecodeRegion(), memoized in #omc_geom_cache against the last region + decoded on this thread. */ static inline void omcCachedDecodeRegion(int irl, int imax, int jmax, int *ix, int *iy, int *iz) { @@ -225,7 +272,8 @@ static inline void omcCachedDecodeRegion(int irl, int imax, int jmax, *iz = omc_geom_cache.iz; } -/* Prefill the memo with a region whose indices the caller already knows */ +/*! Prefill #omc_geom_cache with a region whose indices the caller already + knows. */ static inline void omcStageRegion(int irl, int ix, int iy, int iz) { omc_geom_cache.irl = irl; @@ -234,6 +282,8 @@ static inline void omcStageRegion(int irl, int ix, int iy, int iz) { omc_geom_cache.iz = iz; } +/*! Reciprocals of a direction, memoized in #omc_geom_cache against the last + direction seen on this thread. */ static inline void omcInvDir(double u, double v, double w, double *ui, double *vi, double *wi) { @@ -251,9 +301,8 @@ static inline void omcInvDir(double u, double v, double w, *vi = omc_geom_cache.vi; *wi = omc_geom_cache.wi; } -/******************************************************************************/ -/* Flag set by '--verbose' argument */ +/*! Flag set by the '--verbose' command line argument. */ extern int verbose_flag; -#endif \ No newline at end of file +#endif diff --git a/ucodes/omc_python/ompmc/__init__.py b/ucodes/omc_python/ompmc/__init__.py index 34512fb..81d303a 100644 --- a/ucodes/omc_python/ompmc/__init__.py +++ b/ucodes/omc_python/ompmc/__init__.py @@ -55,6 +55,7 @@ __version__ = _ompmc.__version__ MAX_MEDIA = _ompmc.MAX_MEDIA +"""Maximum number of distinct media a :class:`Geometry` may reference.""" def data_path() -> Path: @@ -63,6 +64,17 @@ def data_path() -> Path: Set ``OMPMC_DATA_PATH`` to override it; otherwise the copy shipped inside the package is used, falling back to the source tree when running from a checkout that was not installed. + + Returns + ------- + pathlib.Path + Directory containing the ``data``, ``pegs4`` and ``spectra`` + subdirectories. + + Raises + ------ + FileNotFoundError + If no such directory can be found and ``OMPMC_DATA_PATH`` is not set. """ override = os.environ.get("OMPMC_DATA_PATH") if override: @@ -103,7 +115,40 @@ def _as_triples(values, name: str) -> np.ndarray: @dataclass class Geometry: - """The voxel phantom: where the boundaries are and what is in each voxel.""" + """The voxel phantom: where the boundaries are and what is in each voxel. + + Parameters + ---------- + x_bounds, y_bounds, z_bounds : array_like + Strictly ascending voxel boundaries along each axis, in cm. ``n + 1`` + values describe ``n`` voxels along that axis. + materials : sequence of str + Medium names, matching entries in the PEGS file. `material` below + indexes into this list, starting from 1; at most :data:`MAX_MEDIA` + are supported. + density : numpy.ndarray + Fortran-ordered ``float64`` cube of mass densities in g/cm^3, shaped + ``(len(x_bounds) - 1, len(y_bounds) - 1, len(z_bounds) - 1)``. + material : numpy.ndarray + Fortran-ordered ``int32`` cube of the same shape, indexing + `materials` from 1; 0 means vacuum. + + Raises + ------ + ValueError + If a bounds vector is not strictly ascending, `materials` is empty + or longer than :data:`MAX_MEDIA`, the cubes are not shaped like the + bounds describe, `density` is negative anywhere, or `material` holds + an index outside ``0 .. len(materials)``. + TypeError + If `density` or `material` do not have the required dtype. + + Notes + ----- + The transport indexes voxels with the first axis varying fastest, so a + C-ordered cube would describe a transposed phantom -- it is rejected + rather than silently copied. Use ``np.asfortranarray(...)``. + """ x_bounds: np.ndarray y_bounds: np.ndarray @@ -164,10 +209,12 @@ def _check_cube(values, dtype, name: str) -> np.ndarray: @property def shape(self) -> tuple[int, int, int]: + """Voxel count along each axis, ``(nx, ny, nz)``.""" return self.density.shape @property def n_voxels(self) -> int: + """Total number of voxels, ``nx * ny * nz``.""" return int(self.density.size) @@ -175,15 +222,30 @@ def n_voxels(self) -> int: class Spectrum: """The energy distribution of the source particles. - Build one with :meth:`from_file`, :meth:`from_histogram` or - :meth:`monoenergetic`. + Build one with :meth:`from_file`, :meth:`from_histogram`, + :meth:`monoenergetic` or :meth:`default`; there is no public constructor. """ _payload: dict @classmethod def from_file(cls, path) -> "Spectrum": - """Read an EGSnrc style ``.spectrum`` file.""" + """Read an EGSnrc style ``.spectrum`` file. + + Parameters + ---------- + path : str or os.PathLike + Path to the spectrum file. + + Returns + ------- + Spectrum + + Raises + ------ + FileNotFoundError + If `path` does not exist. + """ path = Path(path) if not path.is_file(): raise FileNotFoundError(f"No spectrum file at {path}") @@ -192,7 +254,33 @@ def from_file(cls, path) -> "Spectrum": @classmethod def from_histogram(cls, energy, fluence, e_min: float = 0.0, per_mev: bool = False) -> "Spectrum": - """A histogram: ``energy`` holds the upper edge of each bin in MeV.""" + """A histogram spectrum. + + Parameters + ---------- + energy : array_like + Upper edge of each bin, in MeV; finite and strictly ascending. + fluence : array_like + Relative number of particles per bin, same length as `energy`; + non-negative, summing to a positive, finite value. + e_min : float, optional + Lower edge of the first bin, in MeV. Must be less than + ``energy[0]``. + per_mev : bool, optional + If true, `fluence` holds counts per MeV rather than counts per + bin. + + Returns + ------- + Spectrum + + Raises + ------ + ValueError + If the shapes of `energy` and `fluence` disagree, `energy` is + not finite and strictly ascending above `e_min`, or `fluence` + is not non-negative and finite-summing to a positive value. + """ energy = np.ascontiguousarray(energy, dtype=np.float64) fluence = np.ascontiguousarray(fluence, dtype=np.float64) @@ -224,7 +312,22 @@ def from_histogram(cls, energy, fluence, e_min: float = 0.0, @classmethod def monoenergetic(cls, energy: float) -> "Spectrum": - """Every particle starts with the same kinetic energy, in MeV.""" + """Every particle starts with the same kinetic energy. + + Parameters + ---------- + energy : float + Kinetic energy in MeV; positive and finite. + + Returns + ------- + Spectrum + + Raises + ------ + ValueError + If `energy` is not positive and finite. + """ energy = float(energy) if not energy > 0.0 or not np.isfinite(energy): raise ValueError(f"energy is {energy} MeV, it must be positive " @@ -233,7 +336,12 @@ def monoenergetic(cls, energy: float) -> "Spectrum": @classmethod def default(cls) -> "Spectrum": - """The 6 MV bremsstrahlung spectrum shipped with ompMC.""" + """The 6 MV bremsstrahlung spectrum shipped with ompMC. + + Returns + ------- + Spectrum + """ return cls.from_file(data_path() / "spectra" / "mohan6.spectrum") @@ -241,9 +349,24 @@ def default(cls) -> "Spectrum": class BeamletSource: """Beamlet apertures at isocentre, one sparse Dij column each. - ``source`` holds one xyz row per beam, the other three one row per beamlet: - the corner of the aperture rectangle and the two edge vectors spanning it. - ``i_beam`` says which beam each beamlet belongs to, counting from 0. + Parameters + ---------- + i_beam : array_like of int + Index of the beam each beamlet belongs to, counting from 0. + source : array_like + Shape ``(n_beams, 3)``, the xyz source position of each beam. + corner : array_like + Shape ``(n_beamlets, 3)``, the corner of each beamlet's aperture + rectangle. + side1, side2 : array_like + Shape ``(n_beamlets, 3)``, the two edge vectors spanning the + aperture rectangle from `corner`. + + Raises + ------ + ValueError + If `i_beam` is empty, `corner`/`side1`/`side2` do not each have one + row per beamlet, or `i_beam` references a beam outside `source`. """ i_beam: np.ndarray @@ -278,12 +401,27 @@ def __post_init__(self) -> None: @property def n_beamlets(self) -> int: + """Number of beamlets, ``len(i_beam)``.""" return int(self.i_beam.size) @dataclass class CollimatedSource: - """A point source at ``ssd`` behind a rectangular opening on the surface.""" + """A point source at `ssd` behind a rectangular opening on the surface. + + Parameters + ---------- + ssd : float + Distance from the source to the phantom surface, in cm; positive. + x_min, x_max, y_min, y_max : float + Bounds of the rectangular opening on the phantom surface, in cm. + + Raises + ------ + ValueError + If `ssd` is not positive, or the opening has negative width in + either direction. + """ ssd: float x_min: float @@ -300,7 +438,25 @@ def __post_init__(self) -> None: @dataclass class Physics: - """Transport parameters and where the interaction data lives.""" + """Transport parameters and where the interaction data lives. + + Parameters + ---------- + pegs_file, pgs4form_file, data_folder, output_folder : str, os.PathLike or None, optional + Override the corresponding file or directory; each defaults to the + matching path under :func:`data_path`. + global_ecut, global_pcut : float, optional + Global electron and photon transport cut-offs, in MeV. + n_split : int, optional + Photon splitting factor at the source; ``1`` disables splitting. + seeds : tuple of int, optional + The two seeds of the Philox4x32-10 random number generator. + esave : float or None, optional + Electron range-rejection threshold, in MeV; ``None`` disables it. + e_rr, f_rr : float or None, optional + Russian roulette threshold, in MeV, and survival factor. Both must + be set (`f_rr` > 1) to take effect. + """ pegs_file: str | os.PathLike | None = None pgs4form_file: str | os.PathLike | None = None @@ -318,7 +474,12 @@ class Physics: f_rr: float | None = None def input_items(self) -> dict[str, str]: - """The key/value pairs the core library reads its configuration from.""" + """The key/value pairs the core library reads its configuration from. + + Returns + ------- + dict of str to str + """ root = data_path() pegs = self.pegs_file or root / "pegs4" / "700icru.pegs4dat" @@ -379,14 +540,56 @@ def calc_dij( ): """Calculate the dose influence matrix, one sparse column per beamlet. - Returns a ``scipy.sparse.csc_array`` of shape ``(n_voxels, n_beamlets)`` in - Gy per incident particle, or a ``(dose, variance)`` pair when ``variance`` - is true. Voxels below ``rel_dose_threshold`` of the beamlet maximum are - dropped from the column. - - ``progress`` is called with the fraction finished, in [0, 1], once per - batch and once per beamlet. Returning ``False`` stops the calculation and - raises ``KeyboardInterrupt``; Ctrl-C does the same. + Parameters + ---------- + geometry : Geometry + The voxel phantom. + source : BeamletSource + The beamlets to calculate a column for. + spectrum : Spectrum, optional + Source energy spectrum. Defaults to :meth:`Spectrum.default`. + physics : Physics, optional + Transport parameters and data file locations. Defaults to + ``Physics()``. + n_histories : int, optional + Histories simulated per beamlet. + n_batches : int, optional + Statistical batches per beamlet, at least 2, needed for the + uncertainty estimate. + charge : int, optional + Source particle: ``-1`` electrons, ``0`` photons, ``1`` positrons. + rel_dose_threshold : float, optional + Voxels below this fraction of a beamlet's maximum dose are dropped + from its column. In ``[0, 1)``. + gaussian_source : bool, optional + Spread the starting point over the collimator plane instead of a + point source, softening the penumbra. + source_width : float, optional + Standard deviation of the Gaussian source, in cm. Ignored unless + `gaussian_source` is true. + variance : bool, optional + Also return the variance of the mean, per voxel. + progress : callable, optional + Called with the fraction finished, in ``[0, 1]``, once per batch and + once per beamlet. Returning ``False`` stops the calculation. + verbosity : int, optional + Log level passed to the engine. + + Returns + ------- + scipy.sparse.csc_array + Shape ``(geometry.n_voxels, source.n_beamlets)``, dose in Gy per + incident particle. When `variance` is true, a ``(dose, variance)`` + pair of such arrays instead. + + Raises + ------ + ValueError + If `n_batches`, `n_histories`, `charge` or `rel_dose_threshold` are + out of range. + KeyboardInterrupt + If `progress` returned false, or Ctrl-C was pressed, stopping the + calculation before every beamlet was reported. """ from scipy.sparse import csc_array @@ -451,30 +654,79 @@ def calc_forward( """Calculate the dose of a whole weighted set of beamlets, in one cube. This is what ``calc_dij(...) @ weights`` would give, computed directly. - ``weights`` holds one finite, non-negative value per beamlet and is where - the collimation comes in: a blocked beamlet gets 0, an open one its fluence, - a partly transmitting one a fraction of it. Histories go to the beamlets in - proportion to their weight, so a blocked beamlet costs nothing and the run - time no longer grows with the number of beamlets. - - Returns ``(dose, uncertainty)``, both cubes shaped like the phantom. The - dose is in Gy for exactly these weights -- doubling them doubles it -- - unless ``output_dose`` is false, in which case it is the mean deposited - energy. The uncertainty is relative, and 0.9999999 where nothing was - deposited. - - Two things differ from :func:`calc_dij`. ``n_histories`` counts the whole - calculation rather than one beamlet, so multiply it by ``n_beamlets`` to - keep the same statistics. And there is no ``rel_dose_threshold``: it prunes - columns of a sparse matrix, and there is no matrix here -- which is worth - remembering when comparing the two, since it is the ``calc_dij`` result - that is pruned. - - The weights modulate fluence, not spectrum: a beamlet at 0.02 starts 2% of - the particles, with the spectrum unhardened. Attenuation in a collimator, - its scatter and the beam hardening that goes with it are not modelled. - - ``progress`` works as it does for :func:`calc_dij`, called once per batch. + `weights` is where the collimation comes in: a blocked beamlet gets 0, an + open one its fluence, a partly transmitting one a fraction of it. + Histories go to the beamlets in proportion to their weight, so a blocked + beamlet costs nothing and the run time no longer grows with the number of + beamlets. + + Parameters + ---------- + geometry : Geometry + The voxel phantom. + source : BeamletSource + The weighted beamlets. + weights : array_like + One finite, non-negative value per beamlet, summing to a finite, + positive total. Modulates fluence, not spectrum: a beamlet at 0.02 + starts 2% of the particles, with the spectrum unhardened. + Attenuation in a collimator, its scatter and the beam hardening that + goes with it are not modelled. + spectrum : Spectrum, optional + Source energy spectrum. Defaults to :meth:`Spectrum.default`. + physics : Physics, optional + Transport parameters and data file locations. Defaults to + ``Physics()``. + n_histories : int, optional + Histories simulated over the whole calculation -- unlike + :func:`calc_dij`, this does not count per beamlet, so multiply it by + `source.n_beamlets` to keep the same statistics as a `calc_dij` run. + n_batches : int, optional + Statistical batches, at least 2, needed for the uncertainty + estimate. + charge : int, optional + Source particle: ``-1`` electrons, ``0`` photons, ``1`` positrons. + gaussian_source : bool, optional + Spread the starting point over the collimator plane instead of a + point source, softening the penumbra. + source_width : float, optional + Standard deviation of the Gaussian source, in cm. Ignored unless + `gaussian_source` is true. + output_dose : bool, optional + If true, the dose is in Gy for exactly these weights -- doubling + them doubles it. If false, the mean deposited energy is returned + instead. + progress : callable, optional + Called with the fraction finished, in ``[0, 1]``, once per batch. + Returning ``False`` stops the calculation. + verbosity : int, optional + Log level passed to the engine. + + Returns + ------- + dose : numpy.ndarray + Cube shaped like the phantom. + uncertainty : numpy.ndarray + Cube shaped like the phantom, the relative uncertainty of `dose`, + and 0.9999999 where nothing was deposited. + + Raises + ------ + ValueError + If `n_batches`, `n_histories` or `charge` are out of range, `weights` + does not have one entry per beamlet, holds a negative or non-finite + value, or sums to zero or a non-finite value. + KeyboardInterrupt + If `progress` returned false, or Ctrl-C was pressed, before any + result was available -- the batches are averaged, so a run stopped + partway through has no result to return. + + Notes + ----- + There is no ``rel_dose_threshold`` as in :func:`calc_dij`: it prunes + columns of a sparse matrix, and there is no matrix here. Worth + remembering when comparing the two, since it is the `calc_dij` result + that gets pruned. """ _check_run(n_histories, n_batches, charge) @@ -546,12 +798,45 @@ def calc_cube( ): """Calculate the dose everywhere in the phantom from one collimated beam. - Returns ``(dose, uncertainty)``, both cubes shaped like the phantom. The - dose is in Gy per incident fluence unless ``output_dose`` is false, in - which case it is the mean deposited energy. The uncertainty is relative, - and 0.9999999 wherever nothing was deposited. - - ``progress`` works as it does for :func:`calc_dij`, called once per batch. + Parameters + ---------- + geometry : Geometry + The voxel phantom. + source : CollimatedSource + The point source and its collimator opening. + spectrum : Spectrum, optional + Source energy spectrum. Defaults to :meth:`Spectrum.default`. + physics : Physics, optional + Transport parameters and data file locations. Defaults to + ``Physics()``. + n_histories : int, optional + Histories simulated. + n_batches : int, optional + Statistical batches, at least 2, needed for the uncertainty + estimate. + charge : int, optional + Source particle: ``-1`` electrons, ``0`` photons, ``1`` positrons. + output_dose : bool, optional + If true, the dose is in Gy per incident fluence. If false, the mean + deposited energy is returned instead. + progress : callable, optional + Called with the fraction finished, in ``[0, 1]``, once per batch. + Returning ``False`` stops the calculation. + verbosity : int, optional + Log level passed to the engine. + + Returns + ------- + dose : numpy.ndarray + Cube shaped like the phantom. + uncertainty : numpy.ndarray + Cube shaped like the phantom, the relative uncertainty of `dose`, + and 0.9999999 where nothing was deposited. + + Raises + ------ + ValueError + If `n_batches`, `n_histories` or `charge` are out of range. """ _check_run(n_histories, n_batches, charge) From 358e49f41c28b1fadc976f5441ea1b4211452002 Mon Sep 17 00:00:00 2001 From: Niklas Wahl Date: Wed, 5 Aug 2026 13:41:11 +0200 Subject: [PATCH 2/6] add matlab interface documentation --- docs/conf.py | 10 +++ docs/index.md | 13 ++++ docs/matlab-api/index.md | 38 ++++++++++ docs/requirements.txt | 1 + pyproject.toml | 1 + ucodes/omc_matrad/omc_matrad.m | 123 +++++++++++++++++++++++++++++++++ 6 files changed, 186 insertions(+) create mode 100644 docs/matlab-api/index.md create mode 100644 ucodes/omc_matrad/omc_matrad.m diff --git a/docs/conf.py b/docs/conf.py index 3174a52..43a743e 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -50,6 +50,7 @@ def _version_from_cmake() -> str: "sphinx.ext.viewcode", "sphinx.ext.autosectionlabel", "breathe", + "sphinxcontrib.matlab", "sphinx_copybutton", "sphinx_design", ] @@ -172,6 +173,15 @@ def _run_doxygen() -> None: ("py:class", r"sequence"), ] +# -- MATLAB domain ------------------------------------------------------- + +# omc_matrad is a compiled MEX file with no .m source; ucodes/omc_matrad/ +# carries omc_matrad.m as a help-text-only stub next to it (MathWorks' own +# convention for documenting a MEX file -- MATLAB always runs the MEX file +# itself and only reads help/doc from the .m file), which this indexes +# alongside the real recordProgressCallback.m helper. +matlab_src_dir = str(ROOT_DIR / "ucodes" / "omc_matrad") + # -- HTML output ------------------------------------------------------------ html_theme = "furo" diff --git a/docs/index.md b/docs/index.md index bdf8c01..007da54 100644 --- a/docs/index.md +++ b/docs/index.md @@ -28,6 +28,12 @@ Build ompMC and its Python extension. The engine headers a new host embeds ompMC through. ::: +:::{grid-item-card} {octicon}`terminal` MATLAB / Octave +:link: matlab-api/index +:link-type: doc +The `omc_matrad` MEX interface for matRad. +::: + :::{grid-item-card} {octicon}`mark-github` Source :link: https://github.com/e0404/ompMC The repository, issue tracker and full README. @@ -54,3 +60,10 @@ python-api/index c-api/index ``` + +```{toctree} +:hidden: +:caption: MATLAB / Octave + +matlab-api/index +``` diff --git a/docs/matlab-api/index.md b/docs/matlab-api/index.md new file mode 100644 index 0000000..de4522f --- /dev/null +++ b/docs/matlab-api/index.md @@ -0,0 +1,38 @@ +# MATLAB / Octave interface + +`omc_matrad` is a MATLAB/Octave MEX file built from +[ucodes/omc_matrad/omc_matrad.c](https://github.com/e0404/ompMC/blob/master/ucodes/omc_matrad/omc_matrad.c) -- +see {doc}`../getting-started/installation` for how it's built, and note the +`mexLock()` caveat there (a rebuilt MEX file needs a MATLAB restart). + +A MEX file carries no MATLAB source for `help`/`doc` to read, so it ships +with a same-named `.m` file containing only a help comment -- MATLAB always +executes the MEX file and only reads documentation from the `.m` file. That +stub is what the reference below is generated from, so it is also what +`help omc_matrad` prints in MATLAB itself. + +```matlab +addpath('build/bin'); +[dij, dijVar] = omc_matrad(cubeRho, cubeMatIx, mcGeo, mcSrc, mcOpt); +``` + +```{eval-rst} +.. mat:module:: . +``` + +## Reference + +```{eval-rst} +.. mat:autofunction:: omc_matrad +``` + +## Progress callback + +`mcOpt.progressCallback`, if given, replaces the built-in `waitbar`. It must +be a function handle taking a single scalar in `[0, 1]`; the example below is +what the test suite uses to record progress into a variable instead of a +window. + +```{eval-rst} +.. mat:autofunction:: recordProgressCallback +``` diff --git a/docs/requirements.txt b/docs/requirements.txt index 3b70d65..efd4e6e 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -7,6 +7,7 @@ sphinx>=7 furo myst-parser breathe +sphinxcontrib-matlabdomain sphinx-copybutton sphinx-design numpy>=1.22 diff --git a/pyproject.toml b/pyproject.toml index 119e05b..8825048 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -32,6 +32,7 @@ docs = [ "furo", "myst-parser", "breathe", + "sphinxcontrib-matlabdomain", "sphinx-copybutton", "sphinx-design", ] diff --git a/ucodes/omc_matrad/omc_matrad.m b/ucodes/omc_matrad/omc_matrad.m new file mode 100644 index 0000000..6e8a604 --- /dev/null +++ b/ucodes/omc_matrad/omc_matrad.m @@ -0,0 +1,123 @@ +function [dij, dijVar] = omc_matrad(cubeRho, cubeMatIx, mcGeo, mcSrc, mcOpt) %#ok +%OMC_MATRAD Monte Carlo dose calculation for matRad. +% +% [dij, dijVar] = OMC_MATRAD(cubeRho, cubeMatIx, mcGeo, mcSrc, mcOpt) +% +% This is a compiled MEX file (.mexw64/.mexa64/.mexmaca64/...); this .m +% file exists only to document it -- MATLAB reads help text from a +% plain-text .m file of the same name next to a MEX file, but always +% executes the MEX file itself, never this one. See BUILDING.md for how +% it is built. +% +% Returns either a sparse dose-influence matrix (mcOpt.mode = 'dij', +% the default) or the dense dose cube of one weighted field +% (mcOpt.mode = 'forward_beamlet'); see "Mode" below. +% +% version = OMC_MATRAD('version') returns the ompMC version string +% without touching any dose calculation state; OMC_MATRAD('-v') and +% OMC_MATRAD('--version') with no output argument print it instead. +% +% The MEX file locks itself in memory on first use (mexLock) because an +% OpenMP-using MEX file cannot safely be unloaded once a parallel +% region has run. A rebuilt MEX file is therefore only picked up after +% restarting MATLAB; see BUILDING.md. +% +% Inputs: +% +% * cubeRho -- 3D double cube of mass densities, g/cm^3. +% * cubeMatIx -- 3D int32 cube of material indices into mcGeo.material, same size as cubeRho; 0 is not a valid index here (unlike the Python interface, there is no vacuum sentinel). +% * mcGeo -- dose grid, a struct with fields material (cell array of medium names, matching entries in the PEGS file), xBounds, yBounds, zBounds (voxel boundaries along each axis, cm, ascending, one more value than voxels along that axis). +% * mcSrc -- beamlet source, a struct with fields nBixels (number of beamlets), iBeam (index of the beam of each beamlet, counting from 1), xSource/ySource/zSource (source position of each beam), xCorner/yCorner/zCorner (corner of each beamlet's aperture rectangle), xSide1/ySide1/zSide1 and xSide2/ySide2/zSide2 (the two edge vectors of that rectangle), and, only for mcOpt.mode = 'forward_beamlet', bixelWeights (one non-negative weight per beamlet). +% * mcOpt -- run settings, see "Options" below. +% +% Outputs, mode 'dij' (default): +% +% * dij -- sparse double matrix, one column per beamlet, one row per dose-grid voxel, dose in Gy per incident particle. Entries below mcOpt.relDoseThreshold (relative to the column maximum) are dropped. +% * dijVar -- sparse matrix of the same shape and sparsity pattern: variance of the mean, dose units squared. Only computed if requested (nargout >= 2). +% +% Outputs, mode 'forward_beamlet': +% +% * dij -- dense double cube, the size of cubeRho: the dose (or mean deposited energy, see mcOpt.outputDose) of the whole weighted field in mcSrc.bixelWeights, in Gy for exactly those weights. +% * dijVar -- dense cube of the same size: the relative uncertainty of each voxel, 0.9999999 where nothing was deposited -- the convention omc_dosxyz writes into a .3ddose file. Only computed if requested (nargout >= 2). +% +% Example:: +% +% addpath('build/bin'); +% [dij, dijVar] = omc_matrad(cubeRho, cubeMatIx, mcGeo, mcSrc, mcOpt); +% +% Options (fields of mcOpt). nSplit, charge, global_ecut, global_pcut, +% randomSeeds, pegsFile, pgs4formFile, dataFolder and outputFolder are +% required; the rest have defaults. +% +% * nHistories -- histories per beamlet ('dij') or in total ('forward_beamlet'). +% * nBatches -- statistical batches, at least 2. +% * nSplit -- photon splitting factor at the source; 1 disables it. +% * charge -- -1 electrons, 0 photons, +1 positrons. +% * global_ecut -- global electron transport cut-off, MeV. +% * global_pcut -- global photon transport cut-off, MeV. +% * randomSeeds -- [seed1 seed2] for the Philox4x32-10 RNG. +% * pegsFile -- path to the .pegs4dat file. +% * pgs4formFile -- path to the pgs4form.dat file. +% * dataFolder -- path to the cross section data directory. +% * outputFolder -- path ompMC's own diagnostics may be written to. +% * mode -- 'dij' (default) or 'forward_beamlet', see "Mode" below. +% * spectrum -- struct(energy, fluence, eMin, mode), overrides spectrumFile; see "Spectrum" below. +% * spectrumFile -- path to a .spectrum file; default ./spectra/mohan6.spectrum. +% * monoEnergy -- single kinetic energy in MeV, used if neither spectrum nor spectrumFile is given. +% * sourceGeometry -- 'point' (default) or 'gaussian'. +% * sourceGaussianWidth -- standard deviation in cm, 'gaussian' only. +% * relDoseThreshold -- mode 'dij' only: voxels below this fraction of a beamlet's maximum dose are dropped from its column. +% * outputDose -- 1 (default) for dose in Gy, 0 for mean deposited energy. +% * verbose -- 0, 1 or 2; 2 also shows a waitbar unless progressCallback is given. +% * progressCallback -- function handle called with a scalar in [0,1] once per batch and once per finished beamlet, replacing the built-in waitbar; see recordProgressCallback.m for an example. +% * esave, e_rr, f_rr -- variance reduction, see "Variance reduction" below. +% +% Mode: dose-influence matrix or forward dose. 'forward_beamlet' +% computes the dose of a whole weighted field in one go instead of one +% sparse column per beamlet. The collimation is given as +% mcSrc.bixelWeights, a non-negative vector of length mcSrc.nBixels: a +% blocked beamlet gets 0, an open one its fluence, a partly +% transmitting one a fraction of it -- the fluence map matRad already +% optimises, so no new geometry is needed. The result is what dij*w +% would have been, reached directly instead of through the matrix. +% +% Two things change meaning in this mode. First, nHistories counts the +% whole calculation, not one beamlet: switching a 'dij' run over +% unchanged divides the statistics by nBixels, so multiply nHistories +% by nBixels to keep them. Second, relDoseThreshold does nothing: it +% prunes columns of a sparse matrix, and there is no matrix here -- it +% is the 'dij' result that is pruned, so set it to 0 for a +% like-for-like comparison of the two modes. +% +% The weights modulate fluence, not spectrum: a leaf transmitting 2% +% starts 2% of the particles, with the spectrum unhardened. Attenuation +% in the collimator, its scatter and the beam hardening that goes with +% it are not modelled. +% +% Spectrum. The source spectrum is tried in this order: mcOpt.spectrum, +% then mcOpt.spectrumFile, then mcOpt.monoEnergy; whichever loses is +% announced rather than silently dropped. Giving none of them uses +% spectra/mohan6.spectrum. mcOpt.spectrum is a struct with fields +% energy (upper energy of each bin in MeV, strictly ascending), fluence +% (relative number of particles per bin, same length, non-negative), +% eMin (lower energy of the first bin in MeV, optional, default 0) and +% mode (0 for counts per bin, the default, or 1 for counts per MeV). +% Within a bin the energy is sampled uniformly, as it is for a spectrum +% read from file. +% +% Variance reduction: +% +% * nSplit -- uniform photon splitting at the source; > 1 enables it. +% * esave -- electron range rejection: electrons whose residual CSDA range cannot carry them out of the current voxel are terminated below this total energy (MeV). 0 or absent disables it. +% * e_rr, f_rr -- unbiased Russian roulette of newly created electrons below total energy e_rr (MeV), with survival probability 1/f_rr. Both must be set (f_rr > 1) to take effect. +% +% Photon transport uses Woodcock (delta) tracking, so photon steps are +% not stopped at voxel boundaries. +% +% See also RECORDPROGRESSCALLBACK + +% This file is help text only, kept in step with omc_matrad.c by hand: the +% MEX file has no .m source of its own to generate this from. Do not add +% code here -- MATLAB always runs the MEX file, so a body here would never +% execute, and its presence would only make that non-obvious. +end From 804c9b36310f8b2e0dffa5fb218785dc9422abb3 Mon Sep 17 00:00:00 2001 From: Niklas Wahl Date: Wed, 5 Aug 2026 17:05:54 +0200 Subject: [PATCH 3/6] avoid shadowing of mex function with matlab function for mex interface --- .../omc_matrad => docs/_matlab}/omc_matrad.m | 2 +- docs/conf.py | 22 ++++++++++++------ docs/matlab-api/index.md | 23 +++++++++++-------- 3 files changed, 29 insertions(+), 18 deletions(-) rename {ucodes/omc_matrad => docs/_matlab}/omc_matrad.m (99%) diff --git a/ucodes/omc_matrad/omc_matrad.m b/docs/_matlab/omc_matrad.m similarity index 99% rename from ucodes/omc_matrad/omc_matrad.m rename to docs/_matlab/omc_matrad.m index 6e8a604..1dcbcfe 100644 --- a/ucodes/omc_matrad/omc_matrad.m +++ b/docs/_matlab/omc_matrad.m @@ -1,4 +1,4 @@ -function [dij, dijVar] = omc_matrad(cubeRho, cubeMatIx, mcGeo, mcSrc, mcOpt) %#ok +function [dij, dijVar] = omc_matrad(cubeRho, cubeMatIx, mcGeo, mcSrc, mcOpt) %OMC_MATRAD Monte Carlo dose calculation for matRad. % % [dij, dijVar] = OMC_MATRAD(cubeRho, cubeMatIx, mcGeo, mcSrc, mcOpt) diff --git a/docs/conf.py b/docs/conf.py index 43a743e..0d180d9 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -175,17 +175,25 @@ def _run_doxygen() -> None: # -- MATLAB domain ------------------------------------------------------- -# omc_matrad is a compiled MEX file with no .m source; ucodes/omc_matrad/ -# carries omc_matrad.m as a help-text-only stub next to it (MathWorks' own -# convention for documenting a MEX file -- MATLAB always runs the MEX file -# itself and only reads help/doc from the .m file), which this indexes -# alongside the real recordProgressCallback.m helper. -matlab_src_dir = str(ROOT_DIR / "ucodes" / "omc_matrad") +# omc_matrad is a compiled MEX file with no .m source, so docs/_matlab/ +# carries a help-text-only omc_matrad.m stub for this to autodocument. +# It deliberately does NOT live in ucodes/omc_matrad/ alongside the real +# MEX file: build.yml's and CMakeLists.txt's MEX smoke test does +# `addpath('build/bin'); addpath('ucodes/omc_matrad')`, and addpath() +# prepends by default, so a same-named .m file placed there ends up +# ahead of build/bin on the path and permanently shadows the compiled +# MEX -- `exist('omc_matrad', 'file')` stops reporting 3 (MEX-file) and +# every MATLAB/Octave CI job fails. docs/_matlab/ is never added to that +# path, so the stub can only ever be seen by Sphinx. +matlab_src_dir = str(DOCS_DIR / "_matlab") # -- HTML output ------------------------------------------------------------ html_theme = "furo" -html_static_path = ["_static"] +# No custom CSS/JS yet -- an empty _static/ directory is invisible to git +# (it tracks no empty directories) and disappears on a fresh checkout, +# which -W turns into a build failure. Add html_static_path back along +# with the directory once there is an actual asset to put in it. html_title = f"ompMC {version}" html_theme_options = { diff --git a/docs/matlab-api/index.md b/docs/matlab-api/index.md index de4522f..757b2b9 100644 --- a/docs/matlab-api/index.md +++ b/docs/matlab-api/index.md @@ -5,11 +5,13 @@ see {doc}`../getting-started/installation` for how it's built, and note the `mexLock()` caveat there (a rebuilt MEX file needs a MATLAB restart). -A MEX file carries no MATLAB source for `help`/`doc` to read, so it ships -with a same-named `.m` file containing only a help comment -- MATLAB always -executes the MEX file and only reads documentation from the `.m` file. That -stub is what the reference below is generated from, so it is also what -`help omc_matrad` prints in MATLAB itself. +A MEX file carries no MATLAB source for autodoc to read, so the reference +below comes from a same-named `.m` file holding only a help comment -- +the usual way to document a MEX file. It deliberately is **not** shipped +next to the real `omc_matrad` MEX file, though: `addpath` prepends by +default, and a same-named `.m` file anywhere later on the path than +`build/bin` would shadow the compiled MEX file, so `help omc_matrad` in an +actual MATLAB session prints only its default one-liner, not this page. ```matlab addpath('build/bin'); @@ -29,10 +31,11 @@ addpath('build/bin'); ## Progress callback `mcOpt.progressCallback`, if given, replaces the built-in `waitbar`. It must -be a function handle taking a single scalar in `[0, 1]`; the example below is -what the test suite uses to record progress into a variable instead of a -window. +be a function handle taking a single scalar in `[0, 1]`; +[recordProgressCallback.m](https://github.com/e0404/ompMC/blob/master/ucodes/omc_matrad/recordProgressCallback.m) +is what the test suite uses to record progress into a variable instead of a +window: -```{eval-rst} -.. mat:autofunction:: recordProgressCallback +```{literalinclude} ../../ucodes/omc_matrad/recordProgressCallback.m +:language: matlab ``` From 11534717ed913d4d06b0e6de56923cd435660e89 Mon Sep 17 00:00:00 2001 From: Niklas Wahl Date: Wed, 5 Aug 2026 23:59:49 +0200 Subject: [PATCH 4/6] add metadata in pyproject, readme, docs and citation.cff --- CITATION.cff | 42 ++++++++++++++++++++++++++++++++++++++++++ README.md | 2 +- docs/conf.py | 4 ++-- pyproject.toml | 9 ++++++++- 4 files changed, 53 insertions(+), 4 deletions(-) create mode 100644 CITATION.cff diff --git a/CITATION.cff b/CITATION.cff new file mode 100644 index 0000000..77e32c4 --- /dev/null +++ b/CITATION.cff @@ -0,0 +1,42 @@ +# This CITATION.cff file was generated with cffinit. +# Visit https://bit.ly/cffinit to generate yours today! + +cff-version: 1.2.0 +title: ompMC +message: >- + If you use this software, please cite it using the + metadata from this file. +type: software +authors: + - given-names: Edgardo + family-names: Doerner + email: endoerner@gmail.com + affiliation: Pontificia Universidad Católica de Chile + - given-names: Niklas + family-names: Wahl + email: n.wahl@dkfz-heidelberg.de + affiliation: Deutsches Krebsforschungszentrum + orcid: 'https://orcid.org/0000-0002-1451-223X' +repository-code: 'https://github.com/e0404/ompMC' +repository: 'https://github.com/edoerner/ompMC' +keywords: + - radiotherapy + - Monte Carlo +license: GPL-3.0 +preferred-citation: + type: article + title: 'Technical Note: An hybrid parallel implementation for EGSnrc Monte Carlo user codes' + authors: + - given-names: Edgardo + family-names: Doerner + email: endoerner@gmail.com + affiliation: Pontificia Universidad Católica de Chile + - given-names: P. + family-names: Caprile + journal: Medical Physics + volume: 45 + issue: 8 + start: '3969' + end: '3973' + year: 2018 + doi: 10.1002/mp.13033 diff --git a/README.md b/README.md index 23b6e2a..3abcd61 100644 --- a/README.md +++ b/README.md @@ -318,4 +318,4 @@ Windows x64 (MSVC and MinGW), Linux x64, Linux ARM64, macOS x64 and macOS ARM64. ## License GNU General Public License v3.0 — see [LICENSE](LICENSE). -Copyright (C) 2018 Edgardo Doerner (edoerner@fis.puc.cl). +Copyright (C) 2018-2026 Edgardo Doerner and Niklas Wahl. diff --git a/docs/conf.py b/docs/conf.py index 0d180d9..c928d0b 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -25,8 +25,8 @@ # -- Project information ----------------------------------------------------- project = "ompMC" -copyright = "2018, Edgardo Doerner" -author = "Edgardo Doerner" +copyright = "2018-2026, Edgardo Doerner and Niklas Wahl" +author = "Edgardo Doerner and Niklas Wahl" def _version_from_cmake() -> str: diff --git a/pyproject.toml b/pyproject.toml index 8825048..7c45de8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -8,10 +8,17 @@ description = "OpenMP parallel Monte Carlo photon and electron transport in voxe readme = "README.md" requires-python = ">=3.9" license = { file = "LICENSE" } -authors = [{ name = "Edgardo Doerner", email = "edoerner@fis.puc.cl" }] +authors = [ + { name = "Edgardo Doerner", email = "endoerner@gmail.com" }, + { name = "Niklas Wahl", email = "n.wahl@dkfz-heidelberg.de" }, +] +maintainers = [ + { name = "Niklas Wahl", email = "n.wahl@dkfz-heidelberg.de" }, +] classifiers = [ "License :: OSI Approved :: GNU General Public License v3 (GPLv3)", "Programming Language :: C", + "Programming Language :: C++", "Programming Language :: Python :: 3", "Topic :: Scientific/Engineering :: Physics", "Topic :: Scientific/Engineering :: Medical Science Apps.", From fdcb97537b6987d303bc438a32408e2e94eedbab Mon Sep 17 00:00:00 2001 From: Niklas Wahl Date: Thu, 6 Aug 2026 00:00:11 +0200 Subject: [PATCH 5/6] bump version --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 9f9eaee..304cc12 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -21,7 +21,7 @@ cmake_minimum_required(VERSION 3.20) project(ompMC - VERSION 0.1.0 + VERSION 0.2.0 DESCRIPTION "An OpenMP parallel implementation for Monte Carlo particle transport simulations" HOMEPAGE_URL "https://github.com/e0404/ompMC" LANGUAGES C) From b83a9f35496bead3ac0652810e9a1ba241322bb7 Mon Sep 17 00:00:00 2001 From: Niklas Wahl Date: Thu, 6 Aug 2026 00:16:36 +0200 Subject: [PATCH 6/6] update workflows (coverage + dependabot) and badges --- .github/dependabot.yml | 48 +++++++++++++++++++++++++++++++++++++ .github/workflows/build.yml | 6 +++++ README.md | 13 ++++++++++ 3 files changed, 67 insertions(+) diff --git a/.github/dependabot.yml b/.github/dependabot.yml index c71be76..7d9d00a 100644 --- a/.github/dependabot.yml +++ b/.github/dependabot.yml @@ -4,6 +4,14 @@ version: 2 # pull request when a newer release of an action referenced there appears, which # runs the build workflow and therefore gets tested on every runner before it is # merged. +# +# Two more package-ecosystem entries below cover the Python side: the pyproject.toml +# dependencies (numpy/scipy, plus the test/docs extras and the nanobind/scikit-build-core +# build backend) and docs/requirements.txt, which is a separate tree on purpose -- see +# the comment there -- so Dependabot needs to be pointed at it explicitly. +# +# Not covered, because Dependabot has no ecosystem for it: the openlibm FetchContent pin +# (GIT_TAG v0.8.7) in CMakeLists.txt. That one has to be bumped by hand. updates: - package-ecosystem: github-actions @@ -27,3 +35,43 @@ updates: update-types: - minor - patch + + - package-ecosystem: pip + directory: / + schedule: + interval: weekly + day: monday + time: "06:00" + timezone: Europe/Berlin + open-pull-requests-limit: 5 + commit-message: + prefix: deps + labels: + - dependencies + - python + groups: + python-minor: + applies-to: version-updates + update-types: + - minor + - patch + + - package-ecosystem: pip + directory: /docs + schedule: + interval: weekly + day: monday + time: "06:00" + timezone: Europe/Berlin + open-pull-requests-limit: 5 + commit-message: + prefix: docs + labels: + - dependencies + - docs + groups: + docs-minor: + applies-to: version-updates + update-types: + - minor + - patch diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 6c2f510..f7f392f 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -359,3 +359,9 @@ jobs: name: ompmc-coverage path: coverage-html/ if-no-files-found: error + + - name: Upload coverage reports to Codecov + uses: codecov/codecov-action@v5 + with: + token: ${{ secrets.CODECOV_TOKEN }} + files: coverage.info diff --git a/README.md b/README.md index 3abcd61..fbe3273 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,18 @@ # ompMC +[![Build](https://github.com/e0404/ompMC/actions/workflows/build.yml/badge.svg)](https://github.com/e0404/ompMC/actions/workflows/build.yml) +[![Docs](https://github.com/e0404/ompMC/actions/workflows/docs.yml/badge.svg)](https://github.com/e0404/ompMC/actions/workflows/docs.yml) +[![codecov](https://codecov.io/gh/e0404/ompMC/branch/master/graph/badge.svg)](https://codecov.io/gh/e0404/ompMC) +[![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](LICENSE) + +[![C](https://img.shields.io/badge/C-A8B9CC?logo=c&logoColor=white)](src/) +[![Python 3.9+](https://img.shields.io/badge/python-3.9%2B-blue.svg)](pyproject.toml) +[![MATLAB](https://img.shields.io/badge/MATLAB-0076A8?logo=mathworks&logoColor=white)](ucodes/omc_matrad/) + +[![Windows](https://img.shields.io/badge/Windows-0078D6?logo=windows&logoColor=white)](.github/workflows/build.yml) +[![Linux](https://img.shields.io/badge/Linux-FCC624?logo=linux&logoColor=black)](.github/workflows/build.yml) +[![macOS](https://img.shields.io/badge/macOS-000000?logo=apple&logoColor=white)](.github/workflows/build.yml) + > The original repository is **[edoerner/ompMC](https://github.com/edoerner/ompMC)** by Edgardo Doerner. > This repository is a fork under further development, aimed at integration into the > treatment planning toolkits **[matRad](https://github.com/e0404/matRad)** (`e0404/matRad`)