From e88952df9bd5e9165d3a1037b991b21142f39634 Mon Sep 17 00:00:00 2001 From: Tom Body Date: Tue, 21 Apr 2026 17:21:49 -0400 Subject: [PATCH 01/25] Use exact join on datasets --- radas/read_rate_coeffs.py | 87 ++++++++++++++------------------------- 1 file changed, 31 insertions(+), 56 deletions(-) diff --git a/radas/read_rate_coeffs.py b/radas/read_rate_coeffs.py index 3e895fc..18d41b3 100644 --- a/radas/read_rate_coeffs.py +++ b/radas/read_rate_coeffs.py @@ -8,6 +8,8 @@ import xarray as xr import numpy as np +reference_electron_density = Quantity(1.0, ureg.m**-3) +reference_electron_temp = Quantity(1.0, ureg.eV) def read_rate_coeff(data_file_dir, species_name, config): """Builds a rate_dataset combining all of the raw data available for a given species.""" @@ -18,15 +20,7 @@ def read_rate_coeff(data_file_dir, species_name, config): except PackageNotFoundError: radas_version="UNDEFINED" - dataset = xr.Dataset().assign_attrs( - atomic_number=config_for_species["atomic_number"], - species_name=species_name, - git_hash=get_git_revision_short_hash(), - radas_version=radas_version, - created=datetime.date.today().strftime("%Y-%b-%d"), - ) - - dataset = write_global_attributes(dataset, config["globals"]) + rate_coefficients = dict() for dataset_type in config_for_species["data_files"].keys(): reader_key, dataset_config = determine_reader_class_and_config( @@ -40,16 +34,30 @@ def read_rate_coeff(data_file_dir, species_name, config): dataset_type, dataset_config, ) + else: raise NotImplementedError( f"No implementation for reading {reader_key} files." ) - determine_coordinates(dataset, rate_dataset) - dataset[dataset_type] = rate_dataset.rate_coefficient + rate_coefficients[dataset_type] = rate_dataset.rate_coefficient + + dataset = xr.merge([v.rename(k) for k, v in rate_coefficients.items()], join='inner') + dataset["electron_density"] = dataset["dim_electron_density"] * reference_electron_density + dataset["electron_temp"] = dataset["dim_electron_temp"] * reference_electron_temp dataset = align_rates_on_charge_states(dataset) + dataset = dataset.assign_attrs( + atomic_number=config_for_species["atomic_number"], + species_name=species_name, + git_hash=get_git_revision_short_hash(), + radas_version=radas_version, + created=datetime.date.today().strftime("%Y-%b-%d"), + ) + + dataset = write_global_attributes(dataset, config["globals"]) + return dataset @@ -69,26 +77,6 @@ def write_global_attributes(dataset: xr.Dataset, globals: dict) -> xr.Dataset: return dataset -def determine_coordinates(dataset: xr.Dataset, rate_dataset: xr.Dataset): - - for key in [ - "electron_density", - "electron_temp", - "reference_electron_density", - "reference_electron_temp", - ]: - if key not in dataset: - dataset[key] = rate_dataset[key] - - for key in ["electron_density", "electron_temp"]: - np.testing.assert_allclose( - dimensionless_magnitude( - (dataset[key] - rate_dataset[key]) / dataset[f"reference_{key}"] - ), - 0.0, - ) - - def build_adf11_rate_dataset( data_file_dir, species_name, dataset_type, dataset_config ): @@ -107,36 +95,23 @@ def build_adf11_rate_dataset( ) electron_temp = Quantity(10 ** data["DTEVD"][: data["ITMAXD"]], ureg.eV) - # Use logarithmic quantities to define the coordinates, so that we can interpolate over logarithmic quantities. - ds["electron_density"] = xr.DataArray( - electron_density, coords=dict(dim_electron_density=electron_density.magnitude) - ) - ds["electron_temp"] = xr.DataArray( - electron_temp, coords=dict(dim_electron_temp=electron_temp.magnitude) - ) - - ds["reference_electron_density"] = Quantity(1.0, ureg.m**-3) - ds["reference_electron_temp"] = Quantity(1.0, ureg.eV) - - ds["number_of_charge_states"] = data["IZMAX"] - charge_state = np.arange(data["IZMAX"]) - ds["charge_state"] = xr.DataArray( - charge_state, coords=dict(dim_charge_state=charge_state) - ) - coefficient = data["DRCOFD"][: data["IZMAX"], : data["ITMAXD"], : data["IDMAXD"]] if dataset_config["code"] <= 9: coefficient = 10**coefficient input_units = dataset_config["stored_units"] - output_units = dataset_config["desired_units"] - ds["rate_coefficient"] = convert_units( - xr.DataArray( - coefficient, - dims=("dim_charge_state", "dim_electron_temp", "dim_electron_density"), - ).pint.quantify(input_units), - output_units, - ) + + dim_electron_density = dimensionless_magnitude(electron_density / reference_electron_density) + dim_electron_temp = dimensionless_magnitude(electron_temp / reference_electron_temp) + dim_charge_state = np.arange(data["IZMAX"]) + + rate_coefficient = xr.DataArray(coefficient, coords=dict( + dim_charge_state = dim_charge_state, + dim_electron_temp = dim_electron_temp, + dim_electron_density = dim_electron_density, + )).pint.quantify(input_units) + + ds["rate_coefficient"] = convert_units(rate_coefficient, dataset_config["desired_units"]) return ds From 903167420ce9214a9408a7b44fca77fab160786e Mon Sep 17 00:00:00 2001 From: Tom Body Date: Tue, 21 Apr 2026 17:39:07 -0400 Subject: [PATCH 02/25] Add an option to run the full radas command in testing --- .github/workflows/workflow_actions.yml | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/.github/workflows/workflow_actions.yml b/.github/workflows/workflow_actions.yml index 500be92..ae989e4 100644 --- a/.github/workflows/workflow_actions.yml +++ b/.github/workflows/workflow_actions.yml @@ -56,6 +56,13 @@ jobs: radas_config -o ./new_config.yaml radas -s hydrogen -c ./new_config.yaml + - name: Test all species + if: (startsWith(github.ref, 'refs/tags/') || github.event_name == 'workflow_dispatch') && matrix.python-version == '3.12' + env: + MPLBACKEND: Agg + run: | + poetry run radas + # --- JOB 2: BUILD RELEASE ARTIFACTS --- build_release: name: Build Release From 2a3c7e525043d4bf3885402b18671ed9c496e1fb Mon Sep 17 00:00:00 2001 From: Tom Body Date: Tue, 21 Apr 2026 18:38:24 -0400 Subject: [PATCH 03/25] Add option to use inner join for Krypton --- radas/config.yaml | 2 ++ radas/read_rate_coeffs.py | 3 ++- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/radas/config.yaml b/radas/config.yaml index 68734d2..2761c39 100644 --- a/radas/config.yaml +++ b/radas/config.yaml @@ -307,6 +307,8 @@ species: charge_exchange_cross_coupling: 1989 charge_exchange_emission: 1989 mean_ionisation_potential: 1989 + # Ionization/recombination have slightly different electron temp values than line emission. Use inner join to drop missing values, instead of filling in NaNs + join_method: "inner" rubidium: atomic_symbol: "Rb" atomic_number: 37 diff --git a/radas/read_rate_coeffs.py b/radas/read_rate_coeffs.py index 18d41b3..9146857 100644 --- a/radas/read_rate_coeffs.py +++ b/radas/read_rate_coeffs.py @@ -42,7 +42,8 @@ def read_rate_coeff(data_file_dir, species_name, config): rate_coefficients[dataset_type] = rate_dataset.rate_coefficient - dataset = xr.merge([v.rename(k) for k, v in rate_coefficients.items()], join='inner') + dataset = xr.merge([v.rename(k) for k, v in rate_coefficients.items()], + join=config_for_species.get("join_method", "exact")) dataset["electron_density"] = dataset["dim_electron_density"] * reference_electron_density dataset["electron_temp"] = dataset["dim_electron_temp"] * reference_electron_temp From 1fee1f198d3ffdbff1416dfb34c3ddf0531c9077 Mon Sep 17 00:00:00 2001 From: Tom Body Date: Tue, 21 Apr 2026 19:41:06 -0400 Subject: [PATCH 04/25] Tidy up Mavrin plots --- radas/mavrin_reference/compare_to_mavrin.py | 53 +++++++++++++-------- 1 file changed, 32 insertions(+), 21 deletions(-) diff --git a/radas/mavrin_reference/compare_to_mavrin.py b/radas/mavrin_reference/compare_to_mavrin.py index 4f16e88..f7c4974 100644 --- a/radas/mavrin_reference/compare_to_mavrin.py +++ b/radas/mavrin_reference/compare_to_mavrin.py @@ -1,12 +1,11 @@ import xarray as xr -import numpy as np from pathlib import Path import matplotlib.pyplot as plt from .read_mavrin_data import ( read_mavrin_data, compute_Mavrin_polynomial_fit, ) -from ..unit_handling import ureg, magnitude +from ..unit_handling import ureg, magnitude_in_units def compare_radas_to_mavrin(output_dir: Path): @@ -17,7 +16,7 @@ def compare_radas_to_mavrin(output_dir: Path): compare_radas_to_mavrin_per_species(output_dir, species) -def compare_radas_to_mavrin_per_species(output_dir: Path, species: str, max_decades: int = 6): +def compare_radas_to_mavrin_per_species(output_dir: Path, species: str, max_decades: int = 4, show: bool=False): mavrin_data = read_mavrin_data() ds = xr.open_dataset(output_dir / f"{species}.nc").pint.quantify() @@ -46,41 +45,53 @@ def compare_radas_to_mavrin_per_species(output_dir: Path, species: str, max_deca Lz_radas = ds["equilibrium_Lz"].pint.to(ureg.W * ureg.m**3) mean_charge_radas = ds["equilibrium_mean_charge_state"].pint.to(ureg.dimensionless) - fig, axs = plt.subplots(ncols=2) + fig, axs = plt.subplots(ncols=2, nrows=2, sharex="all", sharey="row") for i in range(ds.sizes["dim_ne_tau"]): ne_tau = ds.ne_tau.isel(dim_ne_tau=i).item() - Lz_radas.isel(dim_ne_tau=i).plot(ax=axs[0], label=f"{ne_tau:~P}", color=f"C{i}") + Lz_radas.isel(dim_ne_tau=i).plot(ax=axs[0][0], label=f"{ne_tau:~P}", color=f"C{i}") if Lz_mavrin is not None: - Lz_mavrin.isel(dim_ne_tau=i).plot(ax=axs[0], color=f"C{i}", linestyle="--") + Lz_mavrin.isel(dim_ne_tau=i).plot(ax=axs[0][1], color=f"C{i}") - mean_charge_radas.isel(dim_ne_tau=i).plot(ax=axs[1], color=f"C{i}") + mean_charge_radas.isel(dim_ne_tau=i).plot(ax=axs[1][0], color=f"C{i}") if mean_charge_mavrin is not None: - mean_charge_mavrin.isel(dim_ne_tau=i).plot( - ax=axs[1], color=f"C{i}", linestyle="--" - ) + mean_charge_mavrin.isel(dim_ne_tau=i).plot(ax=axs[1][1], color=f"C{i}") - ds["coronal_Lz"].pint.to(ureg.W * ureg.m**3).plot(ax=axs[0], label="coronal", color="k", linestyle="--") - ds["coronal_mean_charge_state"].pint.to(ureg.dimensionless).plot(ax=axs[1], label="coronal", color="k", linestyle="--") + ds["coronal_Lz"].pint.to(ureg.W * ureg.m**3).plot(ax=axs[0][0], label="coronal", color="k") + ds["coronal_mean_charge_state"].pint.to(ureg.dimensionless).plot(ax=axs[1][0], label="coronal", color="k") - axs[0].legend() - axs[0].set_yscale("log") + axs[0][0].legend() - Lz_radas_mag = magnitude(Lz_radas) - mean_charge_radas_mag = magnitude(mean_charge_radas) + Lz_radas_mag = magnitude_in_units(Lz_radas, ureg.W * ureg.m**3) + Lz_coronal_mag = magnitude_in_units(ds["coronal_Lz"], ureg.W * ureg.m**3) + Lz_min = min(Lz_radas_mag.min(), Lz_coronal_mag.min()) + Lz_max = max(Lz_radas_mag.max(), Lz_coronal_mag.max()) - axs[0].set_ylim(max(np.min(Lz_radas_mag), np.max(Lz_radas_mag) / 10**max_decades) / 2, np.max(Lz_radas_mag) * 2) - axs[1].set_ylim(0, np.max(mean_charge_radas_mag) * 1.2) + for ax in axs.flatten(): + ax.set_title("") + ax.set_ylabel("") + + axs[0][0].set_yscale("log") + axs[0][0].set_ylim(max(Lz_min, Lz_max / 10**max_decades) / 2, Lz_max * 2) + axs[1][0].set_ylim(0, ds.atomic_number * 1.2) + + axs[0][0].set_ylabel("$L_z$ $[W m^3]$") + axs[1][0].set_ylabel("$$") - axs[0].set_title("$L_z$ $[W m^3]$") - axs[1].set_title("$$") + axs[0][0].set_title("radas") + axs[0][1].set_title("Mavrin") for ax in axs.flatten(): ax.set_xscale("log") + ax.set_xlabel("") + + for ax in axs[-1][:].flatten(): ax.set_xlabel("$T_e$ [$eV$]") - ax.set_ylabel("") plt.suptitle(species) + if show: + plt.show() + plt.savefig(output_dir / f"{species}.png") From b3d5afb3a2d6b4aec96a4ccff8c4ec012d62aefe Mon Sep 17 00:00:00 2001 From: Tom Body Date: Tue, 21 Apr 2026 19:43:19 -0400 Subject: [PATCH 05/25] Bump the version to 1!1.0.1 --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 297509a..eb39aa1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,7 +5,7 @@ name = "radas" # used a date-based format (e.g., 2025.6.2), package managers will treat # Semantic Versioning (i.e. 1.0.0) as an "older" release. The epoch ensures that # the new versioning takes precedence over all legacy YYYY.MM.version releases. -version = "1!1.0.0" +version = "1!1.0.1" description = "Plasma radiated power calculated using OpenADAS" authors = ["Commonwealth Fusion Systems"] readme = "README.md" From 1ca504031344d873614c4af38558ac1115e2dbae Mon Sep 17 00:00:00 2001 From: Tom Body Date: Wed, 22 Apr 2026 09:13:27 -0400 Subject: [PATCH 06/25] Remove ipdb from radas command test --- .github/workflows/workflow_actions.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/workflow_actions.yml b/.github/workflows/workflow_actions.yml index ae989e4..2e2a4e6 100644 --- a/.github/workflows/workflow_actions.yml +++ b/.github/workflows/workflow_actions.yml @@ -56,11 +56,12 @@ jobs: radas_config -o ./new_config.yaml radas -s hydrogen -c ./new_config.yaml - - name: Test all species + - name: Test radas command if: (startsWith(github.ref, 'refs/tags/') || github.event_name == 'workflow_dispatch') && matrix.python-version == '3.12' env: MPLBACKEND: Agg run: | + poetry install --without dev poetry run radas # --- JOB 2: BUILD RELEASE ARTIFACTS --- From 9b2c730bc34bc993fb2b5dc2e2d70bca21a0e8fe Mon Sep 17 00:00:00 2001 From: Tom Body Date: Wed, 22 Apr 2026 09:47:22 -0400 Subject: [PATCH 07/25] Prevent ipdb from running without --debug. --- radas/cli.py | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/radas/cli.py b/radas/cli.py index 95a0437..b2ad0be 100644 --- a/radas/cli.py +++ b/radas/cli.py @@ -4,6 +4,7 @@ from pathlib import Path from functools import partial from typing import Optional +import contextlib from .shared import open_yaml_file, default_config_file from .adas_interface.download_adas_datasets import download_species_data @@ -69,15 +70,22 @@ def run_radas_cli( verbose=verbose, debug=debug, ) - try: - from ipdb import launch_ipdb_on_exception - - with launch_ipdb_on_exception(): - run_radas(**kwargs) - except ModuleNotFoundError: + + # nullcontext does nothing. + debug_context = contextlib.nullcontext() + + if debug: + try: + # If --debug and ipdb is installed, switch + # to use the launch_ipdb_on_exception context + from ipdb import launch_ipdb_on_exception + debug_context = launch_ipdb_on_exception() + except ModuleNotFoundError: + pass + + with debug_context: run_radas(**kwargs) - def run_radas( directory: Path, config: Optional[str], From ac33da6b9eb31b9e783541c06ffd63f94b9c8950 Mon Sep 17 00:00:00 2001 From: Tom Body Date: Wed, 22 Apr 2026 09:47:39 -0400 Subject: [PATCH 08/25] Test radas command with 3.12 and upload artifacts --- .github/workflows/workflow_actions.yml | 58 +++++++++++++++++++------- 1 file changed, 44 insertions(+), 14 deletions(-) diff --git a/.github/workflows/workflow_actions.yml b/.github/workflows/workflow_actions.yml index 2e2a4e6..f3be3b4 100644 --- a/.github/workflows/workflow_actions.yml +++ b/.github/workflows/workflow_actions.yml @@ -10,6 +10,7 @@ on: env: FORCE_JAVASCRIPT_ACTIONS_TO_NODE24: true + MPLBACKEND: Agg concurrency: group: ${{ github.workflow }}-${{ github.ref }} @@ -17,7 +18,7 @@ concurrency: jobs: # --- JOB 1: RUN TESTS --- - test: + pytest: name: Test (Python ${{ matrix.python-version }}) runs-on: ubuntu-latest strategy: @@ -41,13 +42,9 @@ jobs: run: poetry install --with dev - name: Run Pytest - env: - MPLBACKEND: Agg run: poetry run pytest - name: Smoke test package installation - env: - MPLBACKEND: Agg run: | poetry build python -m venv test_env @@ -55,19 +52,52 @@ jobs: pip install dist/*.whl radas_config -o ./new_config.yaml radas -s hydrogen -c ./new_config.yaml + + # --- JOB 2: RUN RADAS COMMAND --- + test_radas: + name: Test radas command + runs-on: ubuntu-latest + if: (startsWith(github.ref, 'refs/tags/') || github.event_name == 'workflow_dispatch') + permissions: + contents: write + + steps: + - uses: actions/checkout@v5 + + - name: Install Poetry + run: pipx install "poetry>=2,<3" - - name: Test radas command - if: (startsWith(github.ref, 'refs/tags/') || github.event_name == 'workflow_dispatch') && matrix.python-version == '3.12' + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v6 + with: + python-version: ${{ matrix.python-version }} + cache: 'poetry' + + - name: Install dependencies + run: poetry install + + - name: Run radas command + run: poetry run radas -vvv + + - name: Create and Upload Release Asset + # Only upload if this is a tag (skip for manual workflow_dispatch runs) + if: startsWith(github.ref, 'refs/tags/') env: - MPLBACKEND: Agg + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} run: | - poetry install --without dev - poetry run radas - - # --- JOB 2: BUILD RELEASE ARTIFACTS --- + if [ -d "./radas_dir" ]; then + zip -r radas_dir.zip ./radas_dir + gh release create ${{ github.ref_name }} --generate-notes || true + gh release upload ${{ github.ref_name }} radas_dir.zip --clobber + else + echo "::error::radas_dir not found! Skipping upload." + exit 1 + fi + + # --- JOB 3: BUILD RELEASE ARTIFACTS --- build_release: name: Build Release - needs: test + needs: [pytest, test_radas] if: startsWith(github.ref, 'refs/tags') runs-on: ubuntu-latest steps: @@ -86,7 +116,7 @@ jobs: path: dist/ retention-days: 1 - # --- JOB 3: PUBLISH TO PYPI --- + # --- JOB 4: PUBLISH TO PYPI --- publish: name: Publish to PyPI needs: build_release From 1735a75c8733a00e0d702e72b3d366aa37d1ba43 Mon Sep 17 00:00:00 2001 From: Tom Body Date: Wed, 22 Apr 2026 10:55:53 -0400 Subject: [PATCH 09/25] Record which year data file is from --- radas/adas_interface/download_adas_datasets.py | 2 +- radas/adas_interface/read_adf11_file.py | 6 +++--- radas/read_rate_coeffs.py | 10 +++++++++- 3 files changed, 13 insertions(+), 5 deletions(-) diff --git a/radas/adas_interface/download_adas_datasets.py b/radas/adas_interface/download_adas_datasets.py index a8ec666..c9df09b 100644 --- a/radas/adas_interface/download_adas_datasets.py +++ b/radas/adas_interface/download_adas_datasets.py @@ -32,7 +32,7 @@ def download_species_data( year_key = f"{year}"[-2:] dataset_prefix = dataset_config["prefix"].lower() - output_filename = data_file_dir / f"{species_name}_{dataset_type}.dat" + output_filename = data_file_dir / f"{species_name}_{dataset_type}_{year_key}.dat" query_path = f"{url_base}/download/{reader_class}/{dataset_prefix}{year_key}/{dataset_prefix}{year_key}_{species_key}.dat" if not output_filename.exists(): diff --git a/radas/adas_interface/read_adf11_file.py b/radas/adas_interface/read_adf11_file.py index 3e06fcc..b8075e3 100644 --- a/radas/adas_interface/read_adf11_file.py +++ b/radas/adas_interface/read_adf11_file.py @@ -14,15 +14,15 @@ def load_library(library_name: str, filepath: Path): def read_adf11_file( - data_file_dir, species_name, dataset_type + data_file_dir, species_name, year, dataset_type ) -> dict: """Open and read an ADF11 OpenADAS file. Uses the format specification from https://www.adas.ac.uk/man/appxa-11.pdf """ - - filename = data_file_dir / f"{species_name}_{dataset_type}.dat" + year_key = f"{year}"[-2:] + filename = data_file_dir / f"{species_name}_{dataset_type}_{year_key}.dat" if not filename.exists(): raise FileNotFoundError(f"{filename} does not exist.") diff --git a/radas/read_rate_coeffs.py b/radas/read_rate_coeffs.py index 9146857..994acde 100644 --- a/radas/read_rate_coeffs.py +++ b/radas/read_rate_coeffs.py @@ -22,15 +22,23 @@ def read_rate_coeff(data_file_dir, species_name, config): rate_coefficients = dict() - for dataset_type in config_for_species["data_files"].keys(): + for dataset_type, file_to_download in config_for_species["data_files"].items(): reader_key, dataset_config = determine_reader_class_and_config( config["data_file_config"], dataset_type ) + if isinstance(file_to_download, int): + year = file_to_download + elif isinstance(file_to_download, list) and len(file_to_download) == 2: + year = file_to_download[1] + else: + raise NotImplementedError(f"Could not process entry: {file_to_download} for {species_name} {dataset_type}") + if reader_key == "adf11": rate_dataset = build_adf11_rate_dataset( data_file_dir, species_name, + year, dataset_type, dataset_config, ) From 82b25fa9fb431e9d26db2a09fab0f9e645d747f5 Mon Sep 17 00:00:00 2001 From: Tom Body Date: Wed, 22 Apr 2026 10:56:12 -0400 Subject: [PATCH 10/25] Improve error handling for AlignmentErrors --- radas/read_rate_coeffs.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/radas/read_rate_coeffs.py b/radas/read_rate_coeffs.py index 994acde..9828800 100644 --- a/radas/read_rate_coeffs.py +++ b/radas/read_rate_coeffs.py @@ -50,8 +50,16 @@ def read_rate_coeff(data_file_dir, species_name, config): rate_coefficients[dataset_type] = rate_dataset.rate_coefficient - dataset = xr.merge([v.rename(k) for k, v in rate_coefficients.items()], - join=config_for_species.get("join_method", "exact")) + join_method = config_for_species.get("join_method", "exact") + try: + dataset = xr.merge([v.rename(k) for k, v in rate_coefficients.items()], + join=join_method) + except xr.AlignmentError as e: + raise xr.AlignmentError(f"Alignment failed for {species_name} with join={join_method}. Error was {e}") + + if dataset.sizes["dim_electron_density"] <= 2 or dataset.sizes["dim_electron_temp"] <= 2: + raise xr.AlignmentError(f"Alignment resulted in zero-length axes for {species_name} with join={join_method}. Resulting dataset sizes {dataset.sizes}") + dataset["electron_density"] = dataset["dim_electron_density"] * reference_electron_density dataset["electron_temp"] = dataset["dim_electron_temp"] * reference_electron_temp @@ -87,11 +95,11 @@ def write_global_attributes(dataset: xr.Dataset, globals: dict) -> xr.Dataset: def build_adf11_rate_dataset( - data_file_dir, species_name, dataset_type, dataset_config + data_file_dir, species_name, year, dataset_type, dataset_config ): from .adas_interface.read_adf11_file import read_adf11_file - data = read_adf11_file(data_file_dir, species_name, dataset_type) + data = read_adf11_file(data_file_dir, species_name, year, dataset_type) ds = xr.Dataset() From cd7877fbe61f60055ba179ba5dca55a42132e0f6 Mon Sep 17 00:00:00 2001 From: Tom Body Date: Wed, 22 Apr 2026 10:56:59 -0400 Subject: [PATCH 11/25] Ensure all data files come from the same year to resolve alignment error --- radas/config.yaml | 40 ++++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/radas/config.yaml b/radas/config.yaml index 2761c39..378da6f 100644 --- a/radas/config.yaml +++ b/radas/config.yaml @@ -137,24 +137,24 @@ species: atomic_symbol: "Li" atomic_number: 3 data_files: - effective_recombination: 1996 - effective_ionisation: 1996 - line_emission_from_excitation: 1996 - recombination_and_bremsstrahlung: 1996 + effective_recombination: 1989 + effective_ionisation: 1989 + line_emission_from_excitation: 1989 + recombination_and_bremsstrahlung: 1989 charge_exchange_cross_coupling: 1989 charge_exchange_emission: 1989 - mean_ionisation_potential: 1996 + mean_ionisation_potential: 1989 beryllium: atomic_symbol: "Be" atomic_number: 4 data_files: - effective_recombination: 1996 - effective_ionisation: 1996 - line_emission_from_excitation: 1996 - recombination_and_bremsstrahlung: 1996 + effective_recombination: 1989 + effective_ionisation: 1989 + line_emission_from_excitation: 1989 + recombination_and_bremsstrahlung: 1989 charge_exchange_cross_coupling: 1989 charge_exchange_emission: 1989 - mean_ionisation_potential: 1996 + mean_ionisation_potential: 1989 boron: atomic_symbol: "B" atomic_number: 5 @@ -192,13 +192,13 @@ species: atomic_symbol: "O" atomic_number: 8 data_files: - effective_recombination: 1996 - effective_ionisation: 1996 - line_emission_from_excitation: 1996 - recombination_and_bremsstrahlung: 1996 + effective_recombination: 1989 + effective_ionisation: 1989 + line_emission_from_excitation: 1989 + recombination_and_bremsstrahlung: 1989 charge_exchange_cross_coupling: 1989 charge_exchange_emission: 1989 - mean_ionisation_potential: 1996 + mean_ionisation_potential: 1989 fluorine: atomic_symbol: "F" atomic_number: 9 @@ -206,13 +206,13 @@ species: atomic_symbol: "Ne" atomic_number: 10 data_files: - effective_recombination: 1996 - effective_ionisation: 1996 - line_emission_from_excitation: 1996 - recombination_and_bremsstrahlung: 1996 + effective_recombination: 1989 + effective_ionisation: 1989 + line_emission_from_excitation: 1989 + recombination_and_bremsstrahlung: 1989 charge_exchange_cross_coupling: 1989 charge_exchange_emission: 1989 - mean_ionisation_potential: 1996 + mean_ionisation_potential: 1989 sodium: atomic_symbol: "Na" atomic_number: 11 From d955325450fc19e4be2709c6da505014d5553e6a Mon Sep 17 00:00:00 2001 From: Tom Body Date: Wed, 22 Apr 2026 11:02:03 -0400 Subject: [PATCH 12/25] Fix missing matrix variable in Github actions --- .github/workflows/workflow_actions.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/workflow_actions.yml b/.github/workflows/workflow_actions.yml index f3be3b4..18abdc0 100644 --- a/.github/workflows/workflow_actions.yml +++ b/.github/workflows/workflow_actions.yml @@ -67,10 +67,10 @@ jobs: - name: Install Poetry run: pipx install "poetry>=2,<3" - - name: Set up Python ${{ matrix.python-version }} + - name: Set up Python uses: actions/setup-python@v6 with: - python-version: ${{ matrix.python-version }} + python-version: '3.12' cache: 'poetry' - name: Install dependencies From 9d2cd134b3aed07366950ec1abf100efa514ef81 Mon Sep 17 00:00:00 2001 From: Tom Body Date: Wed, 22 Apr 2026 13:58:20 -0400 Subject: [PATCH 13/25] Use nearest-neighbour if extrapolation needed --- radas/interpolate_rates.py | 54 +++++++------------------------------- 1 file changed, 9 insertions(+), 45 deletions(-) diff --git a/radas/interpolate_rates.py b/radas/interpolate_rates.py index 73b0f57..6b45a7c 100644 --- a/radas/interpolate_rates.py +++ b/radas/interpolate_rates.py @@ -8,6 +8,8 @@ def interpolate_array(array: xr.DataArray, new_electron_density: NDArray[np.floa """Interpolate array onto new values for the electron density and electron temp. The interpolation is performed for logarithmic values. + + Nearest-neighbour extrapolation is used to fill any off-grid points. """ units = array.pint.units array = array.pint.dequantify().squeeze() @@ -27,52 +29,14 @@ def interpolate_array(array: xr.DataArray, new_electron_density: NDArray[np.floa x_interp = np.log10(new_electron_density) y_interp = np.log10(new_electron_temp) - z_interp = np.power(10, RectBivariateSpline(x, y, z)(x_interp, y_interp, grid=True).T) + + # Clip the interpolation to exclude off-grid points. This is equivalent to using nearest-neighbour + # extrapolation for off-grid values. + x_interp_clipped = np.clip(x_interp, x.min().item(), x.max().item()) + y_interp_clipped = np.clip(y_interp, y.min().item(), y.max().item()) + + z_interp = np.power(10, RectBivariateSpline(x, y, z)(x_interp_clipped, y_interp_clipped, grid=True).T) return xr.DataArray(z_interp, coords=dict(dim_electron_temp=new_electron_temp, dim_electron_density=new_electron_density) ) * units - -def interpolate_dataset(dataset: xr.Dataset, electron_density_resolution: int, electron_temp_resolution: int) -> xr.Dataset: - """Interpolate all rate coefficients in a dataset.""" - new_electron_density = np.logspace( - np.log10(dataset["dim_electron_density"].min().item()), - np.log10(dataset["dim_electron_density"].max().item()), - num = electron_density_resolution - ) - - new_electron_temp = np.logspace( - np.log10(dataset["dim_electron_temp"].min().item()), - np.log10(dataset["dim_electron_temp"].max().item()), - num = electron_temp_resolution - ) - - new_dataset = xr.Dataset().assign_attrs(dataset.attrs) - - for key, array in dataset.items(): - if key in [ - "electron_density", - "electron_temp", - ]: - # Don't copy in the coordinate arrays which we'll interpolate - continue - elif key in [ - "ne_tau" - ]: - # Directly copy in the coordinate arrays which we'll leave unchanged - new_dataset[key] = array - elif array.ndim == 0: - # Directly copy in scalar arrays - new_dataset[key] = array - elif (("dim_electron_density" in array.coords) - and ("dim_electron_temp" in array.coords) - and ("dim_charge_state" in array.coords)): - # For each charge state, interpolate the rate coefficient - new_dataset[key] = array.groupby("dim_charge_state").map(interpolate_array, args=(new_electron_density, new_electron_temp)) - else: - raise NotImplementedError(f"Could not process array '{key}' with coords {array.coords}") - - new_dataset["electron_density"] = xr.DataArray(new_electron_density, dims="dim_electron_density") * dataset["electron_density"].pint.units - new_dataset["electron_temp"] = xr.DataArray(new_electron_temp, dims="dim_electron_temp") * dataset["electron_temp"].pint.units - - return new_dataset \ No newline at end of file From 7a2c2f0b097fff0a8f7f9088b255a65bae5d514b Mon Sep 17 00:00:00 2001 From: Tom Body Date: Wed, 22 Apr 2026 13:58:47 -0400 Subject: [PATCH 14/25] Remove optional interpolation of rates. Now required to align grids --- radas/cli.py | 14 -------------- radas/config.yaml | 11 +++++------ 2 files changed, 5 insertions(+), 20 deletions(-) diff --git a/radas/cli.py b/radas/cli.py index b2ad0be..da831b9 100644 --- a/radas/cli.py +++ b/radas/cli.py @@ -15,7 +15,6 @@ from .time_evolution import calculate_time_evolution from .unit_handling import convert_units, ureg from .mavrin_reference import compare_radas_to_mavrin -from .interpolate_rates import interpolate_dataset @click.command() @@ -137,19 +136,6 @@ def run_radas( data_file_dir, species_name, configuration ) - if ("electron_density_resolution" in configuration["globals"]) or ("electron_temp_resolution") in configuration["globals"]: - if verbose: - print("Interpolating rate coefficients") - - new_datasets = dict() - for species_name, dataset in datasets.items(): - electron_density_resolution = configuration["globals"].get("electron_density_resolution", dataset.sizes["dim_electron_density"]) - electron_temp_resolution = configuration["globals"].get("electron_temp_resolution", dataset.sizes["dim_electron_temp"]) - new_datasets[species_name] = interpolate_dataset(dataset, - electron_density_resolution = electron_density_resolution, - electron_temp_resolution = electron_temp_resolution) - datasets = new_datasets - output_dir.mkdir(exist_ok=True, parents=True) if not debug: with mp.Pool() as pool: diff --git a/radas/config.yaml b/radas/config.yaml index 378da6f..77dda1f 100644 --- a/radas/config.yaml +++ b/radas/config.yaml @@ -15,7 +15,6 @@ globals: # Number of log-spaced points to use for electron density and electron temp. # Interpolation is used to map from raw data to points. - # To disable interpolation, do not provide these inputs. electron_density_resolution: 50 electron_temp_resolution: 100 @@ -206,13 +205,13 @@ species: atomic_symbol: "Ne" atomic_number: 10 data_files: - effective_recombination: 1989 - effective_ionisation: 1989 - line_emission_from_excitation: 1989 - recombination_and_bremsstrahlung: 1989 + effective_recombination: 1996 + effective_ionisation: 1996 + line_emission_from_excitation: 1996 + recombination_and_bremsstrahlung: 1996 charge_exchange_cross_coupling: 1989 charge_exchange_emission: 1989 - mean_ionisation_potential: 1989 + mean_ionisation_potential: 1996 sodium: atomic_symbol: "Na" atomic_number: 11 From 54290f3bd5da59477d9e98579ceb4c8cb3534699 Mon Sep 17 00:00:00 2001 From: Tom Body Date: Wed, 22 Apr 2026 13:59:15 -0400 Subject: [PATCH 15/25] Implement interpolation to align data grids --- radas/read_rate_coeffs.py | 88 +++++++++++++++++++++++++++------------ 1 file changed, 61 insertions(+), 27 deletions(-) diff --git a/radas/read_rate_coeffs.py b/radas/read_rate_coeffs.py index 9828800..6bd8d67 100644 --- a/radas/read_rate_coeffs.py +++ b/radas/read_rate_coeffs.py @@ -7,38 +7,69 @@ import datetime import xarray as xr import numpy as np +from .interpolate_rates import interpolate_array reference_electron_density = Quantity(1.0, ureg.m**-3) reference_electron_temp = Quantity(1.0, ureg.eV) def read_rate_coeff(data_file_dir, species_name, config): """Builds a rate_dataset combining all of the raw data available for a given species.""" - config_for_species = config["species"][species_name] try: radas_version=version("radas") except PackageNotFoundError: radas_version="UNDEFINED" + rate_coefficients = build_sorted_dictionary_of_rate_coefficients(config, species_name, data_file_dir) + rate_coefficients = interpolate_rates_onto_matching_grids(config, rate_coefficients) + + try: + dataset = xr.merge([v.rename(k) for k, v in rate_coefficients.items()], join="exact") + except xr.AlignmentError as e: + raise xr.AlignmentError(f"Alignment failed for {species_name} with join='exact'. Error was {e}") + + if dataset.sizes["dim_electron_density"] <= 2 or dataset.sizes["dim_electron_temp"] <= 2: + raise xr.AlignmentError(f"Alignment resulted in dataset sizes {dataset.sizes} for {species_name}.") + + dataset["electron_density"] = dataset["dim_electron_density"] * reference_electron_density + dataset["electron_temp"] = dataset["dim_electron_temp"] * reference_electron_temp + + dataset = align_rates_on_charge_states(dataset) + + dataset = dataset.assign_attrs( + atomic_number=config["species"][species_name]["atomic_number"], + species_name=species_name, + git_hash=get_git_revision_short_hash(), + radas_version=radas_version, + created=datetime.date.today().strftime("%Y-%b-%d"), + ) + + dataset = write_global_attributes(dataset, config["globals"]) + + return dataset + +def build_sorted_dictionary_of_rate_coefficients(config, species_name, data_file_dir): + """Parses configuration data to build a collection of rate coefficients sorted by year (most-recent first).""" rate_coefficients = dict() + years = dict() - for dataset_type, file_to_download in config_for_species["data_files"].items(): + for dataset_type, file_to_read in config["species"][species_name]["data_files"].items(): reader_key, dataset_config = determine_reader_class_and_config( config["data_file_config"], dataset_type ) - if isinstance(file_to_download, int): - year = file_to_download - elif isinstance(file_to_download, list) and len(file_to_download) == 2: - year = file_to_download[1] + if isinstance(file_to_read, int): + years[dataset_type] = file_to_read + elif isinstance(file_to_read, list) and len(file_to_read) == 2: + years[dataset_type] = file_to_read[1] else: - raise NotImplementedError(f"Could not process entry: {file_to_download} for {species_name} {dataset_type}") + raise NotImplementedError(f"Could not process entry: {file_to_read} for {species_name} {dataset_type}") if reader_key == "adf11": rate_dataset = build_adf11_rate_dataset( data_file_dir, species_name, - year, + years[dataset_type], dataset_type, dataset_config, ) @@ -50,32 +81,35 @@ def read_rate_coeff(data_file_dir, species_name, config): rate_coefficients[dataset_type] = rate_dataset.rate_coefficient - join_method = config_for_species.get("join_method", "exact") - try: - dataset = xr.merge([v.rename(k) for k, v in rate_coefficients.items()], - join=join_method) - except xr.AlignmentError as e: - raise xr.AlignmentError(f"Alignment failed for {species_name} with join={join_method}. Error was {e}") + # Sort the datasets so that the most recent data comes first + sorted_by_year = dict(sorted(years.items(), key=lambda item: item[1], reverse=True)) - if dataset.sizes["dim_electron_density"] <= 2 or dataset.sizes["dim_electron_temp"] <= 2: - raise xr.AlignmentError(f"Alignment resulted in zero-length axes for {species_name} with join={join_method}. Resulting dataset sizes {dataset.sizes}") + return {k: rate_coefficients[k] for k in sorted_by_year.keys()} - dataset["electron_density"] = dataset["dim_electron_density"] * reference_electron_density - dataset["electron_temp"] = dataset["dim_electron_temp"] * reference_electron_temp +def interpolate_rates_onto_matching_grids(config, rate_coefficients): - dataset = align_rates_on_charge_states(dataset) + # Since the rate coefficients are sorted by year, taking the first value + # gives the most recent data. + most_recent_rate_coeff = list(rate_coefficients.values())[0] - dataset = dataset.assign_attrs( - atomic_number=config_for_species["atomic_number"], - species_name=species_name, - git_hash=get_git_revision_short_hash(), - radas_version=radas_version, - created=datetime.date.today().strftime("%Y-%b-%d"), + new_electron_density = np.logspace( + np.log10(most_recent_rate_coeff["dim_electron_density"].min().item()), + np.log10(most_recent_rate_coeff["dim_electron_density"].max().item()), + num = config["globals"]["electron_density_resolution"] ) - dataset = write_global_attributes(dataset, config["globals"]) + new_electron_temp = np.logspace( + np.log10(most_recent_rate_coeff["dim_electron_temp"].min().item()), + np.log10(most_recent_rate_coeff["dim_electron_temp"].max().item()), + num = config["globals"]["electron_temp_resolution"] + ) - return dataset + interpolated_rate_coefficients = dict() + + for key, value in rate_coefficients.items(): + interpolated_rate_coefficients[key] = value.groupby("dim_charge_state").map(interpolate_array, args=(new_electron_density, new_electron_temp)) + + return interpolated_rate_coefficients def write_global_attributes(dataset: xr.Dataset, globals: dict) -> xr.Dataset: From 4c5e21b032010481fe53e132c7494b6157ee8784 Mon Sep 17 00:00:00 2001 From: Tom Body Date: Wed, 22 Apr 2026 14:00:06 -0400 Subject: [PATCH 16/25] Revert config to use latest available data --- radas/config.yaml | 34 ++++++++++++++++------------------ 1 file changed, 16 insertions(+), 18 deletions(-) diff --git a/radas/config.yaml b/radas/config.yaml index 77dda1f..4ccc42f 100644 --- a/radas/config.yaml +++ b/radas/config.yaml @@ -136,24 +136,24 @@ species: atomic_symbol: "Li" atomic_number: 3 data_files: - effective_recombination: 1989 - effective_ionisation: 1989 - line_emission_from_excitation: 1989 - recombination_and_bremsstrahlung: 1989 + effective_recombination: 1996 + effective_ionisation: 1996 + line_emission_from_excitation: 1996 + recombination_and_bremsstrahlung: 1996 charge_exchange_cross_coupling: 1989 charge_exchange_emission: 1989 - mean_ionisation_potential: 1989 + mean_ionisation_potential: 1996 beryllium: atomic_symbol: "Be" atomic_number: 4 data_files: - effective_recombination: 1989 - effective_ionisation: 1989 - line_emission_from_excitation: 1989 - recombination_and_bremsstrahlung: 1989 + effective_recombination: 1996 + effective_ionisation: 1996 + line_emission_from_excitation: 1996 + recombination_and_bremsstrahlung: 1996 charge_exchange_cross_coupling: 1989 charge_exchange_emission: 1989 - mean_ionisation_potential: 1989 + mean_ionisation_potential: 1996 boron: atomic_symbol: "B" atomic_number: 5 @@ -191,13 +191,13 @@ species: atomic_symbol: "O" atomic_number: 8 data_files: - effective_recombination: 1989 - effective_ionisation: 1989 - line_emission_from_excitation: 1989 - recombination_and_bremsstrahlung: 1989 + effective_recombination: 1996 + effective_ionisation: 1996 + line_emission_from_excitation: 1996 + recombination_and_bremsstrahlung: 1996 charge_exchange_cross_coupling: 1989 charge_exchange_emission: 1989 - mean_ionisation_potential: 1989 + mean_ionisation_potential: 1996 fluorine: atomic_symbol: "F" atomic_number: 9 @@ -306,8 +306,6 @@ species: charge_exchange_cross_coupling: 1989 charge_exchange_emission: 1989 mean_ionisation_potential: 1989 - # Ionization/recombination have slightly different electron temp values than line emission. Use inner join to drop missing values, instead of filling in NaNs - join_method: "inner" rubidium: atomic_symbol: "Rb" atomic_number: 37 @@ -460,4 +458,4 @@ species: atomic_number: 81 lead: atomic_symbol: "Pb" - atomic_number: 82 + atomic_number: 82 \ No newline at end of file From 2b35e7d84da465796707bd96be55d43d624d936f Mon Sep 17 00:00:00 2001 From: Tom Body Date: Wed, 22 Apr 2026 14:06:37 -0400 Subject: [PATCH 17/25] Run heavier elements first in parallel --- radas/cli.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/radas/cli.py b/radas/cli.py index da831b9..aec68e9 100644 --- a/radas/cli.py +++ b/radas/cli.py @@ -144,11 +144,15 @@ def run_radas( species_name: datasets[species_name] for species_name in species } + # Sort by atomic number, to process heavier elements first since they take longer + # N.b. there will be a race condition, so it might not appear these start first + sorted_datasets = dict(sorted(datasets.items(), key=lambda item: item[1].atomic_number, reverse=True)) + pool.map( partial( run_radas_computation, output_dir=output_dir, verbose=verbose ), - [(ds) for ds in datasets.values()], + [(ds) for ds in sorted_datasets.values()], ) else: for ds in datasets.values(): From 8cac4367a435f264a06e6ff0bfec675f5b27455d Mon Sep 17 00:00:00 2001 From: Tom Body Date: Wed, 22 Apr 2026 14:15:21 -0400 Subject: [PATCH 18/25] Improve documentation of rate reading and interpolation --- radas/cli.py | 3 + radas/interpolate_rates.py | 45 +++++++------ radas/read_rate_coeffs.py | 128 ++++++++++++++++++------------------- 3 files changed, 91 insertions(+), 85 deletions(-) diff --git a/radas/cli.py b/radas/cli.py index aec68e9..f501043 100644 --- a/radas/cli.py +++ b/radas/cli.py @@ -197,6 +197,9 @@ def run_radas_computation(dataset: xr.Dataset, output_dir: Path, verbose: int): output_dir.mkdir(exist_ok=True) dataset.pint.dequantify().to_netcdf(output_dir / f"{dataset.species_name}.nc") + if verbose: + print(f"Finished computation for {dataset.species_name}") + @click.command() @click.option( diff --git a/radas/interpolate_rates.py b/radas/interpolate_rates.py index 6b45a7c..2af7c34 100644 --- a/radas/interpolate_rates.py +++ b/radas/interpolate_rates.py @@ -1,42 +1,51 @@ -"""Routines to interpolate a dataset of rate coefficients to higher resolution.""" +"""Routines for log-log interpolation of rate coefficients with boundary clipping.""" import xarray as xr import numpy as np from scipy.interpolate import RectBivariateSpline from numpy.typing import NDArray -def interpolate_array(array: xr.DataArray, new_electron_density: NDArray[np.floating], new_electron_temp: NDArray[np.floating]) -> xr.DataArray: - """Interpolate array onto new values for the electron density and electron temp. +def interpolate_array( + array: xr.DataArray, + new_electron_density: NDArray[np.floating], + new_electron_temp: NDArray[np.floating] +) -> xr.DataArray: + """ + Interpolate rate coefficients onto a new density/temperature grid in log-log space. - The interpolation is performed for logarithmic values. - - Nearest-neighbour extrapolation is used to fill any off-grid points. + Uses nearest-neighbor extrapolation by clipping out-of-bounds coordinates to + the original grid edges. """ units = array.pint.units array = array.pint.dequantify().squeeze() + # Handle zero-value edge cases (log of zero is undefined) if np.allclose(array, 0.0, atol=0.0, rtol=1e-6): - # If all values of the array are zero, return a zero array. - return xr.DataArray(np.zeros((np.size(new_electron_temp), np.size(new_electron_density))), + return xr.DataArray( + np.zeros((np.size(new_electron_temp), np.size(new_electron_density))), coords=dict(dim_electron_temp=new_electron_temp, dim_electron_density=new_electron_density) ) * units - elif np.any(array <= 0.0): - # If only some of the values of the array are zero, raise an error. - raise NotImplementedError("Cannot handle zero-valued entries in non-zero rate coefficients.") + if np.any(array <= 0.0): + raise NotImplementedError("Cannot log-interpolate rate coefficients containing zeros.") + + # Prepare original grid and data in log10 space x = np.log10(array.dim_electron_density) y = np.log10(array.dim_electron_temp) z = np.log10(array.transpose("dim_electron_density", "dim_electron_temp").pint.magnitude) + # Transform target coordinates to log10 x_interp = np.log10(new_electron_density) y_interp = np.log10(new_electron_temp) - # Clip the interpolation to exclude off-grid points. This is equivalent to using nearest-neighbour - # extrapolation for off-grid values. - x_interp_clipped = np.clip(x_interp, x.min().item(), x.max().item()) - y_interp_clipped = np.clip(y_interp, y.min().item(), y.max().item()) + # Force nearest-neighbor extrapolation by clipping points to the grid domain + x_clipped = np.clip(x_interp, x.min().item(), x.max().item()) + y_clipped = np.clip(y_interp, y.min().item(), y.max().item()) - z_interp = np.power(10, RectBivariateSpline(x, y, z)(x_interp_clipped, y_interp_clipped, grid=True).T) + # Perform spline interpolation and revert from log space + z_interp_log = RectBivariateSpline(x, y, z)(x_clipped, y_clipped, grid=True) + z_interp = np.power(10, z_interp_log.T) - return xr.DataArray(z_interp, + return xr.DataArray( + z_interp, coords=dict(dim_electron_temp=new_electron_temp, dim_electron_density=new_electron_density) - ) * units + ) * units \ No newline at end of file diff --git a/radas/read_rate_coeffs.py b/radas/read_rate_coeffs.py index 6bd8d67..82267bd 100644 --- a/radas/read_rate_coeffs.py +++ b/radas/read_rate_coeffs.py @@ -9,33 +9,43 @@ import numpy as np from .interpolate_rates import interpolate_array +# Reference units for non-dimensionalizing coordinates reference_electron_density = Quantity(1.0, ureg.m**-3) reference_electron_temp = Quantity(1.0, ureg.eV) def read_rate_coeff(data_file_dir, species_name, config): - """Builds a rate_dataset combining all of the raw data available for a given species.""" - + """ + Main pipeline to assemble an atomic rate dataset for a specific species. + + Reads raw ADAS files, standardizes their grids, aligns charge states, + and attaches metadata. + """ try: - radas_version=version("radas") + radas_version = version("radas") except PackageNotFoundError: - radas_version="UNDEFINED" + radas_version = "UNDEFINED" + # 1. Collect and sort data by year rate_coefficients = build_sorted_dictionary_of_rate_coefficients(config, species_name, data_file_dir) + + # 2. Resample all datasets to a common resolution rate_coefficients = interpolate_rates_onto_matching_grids(config, rate_coefficients) + # 3. Merge individual datasets (e.g., recombination, ionization) into one try: dataset = xr.merge([v.rename(k) for k, v in rate_coefficients.items()], join="exact") except xr.AlignmentError as e: - raise xr.AlignmentError(f"Alignment failed for {species_name} with join='exact'. Error was {e}") + raise xr.AlignmentError(f"Alignment failed for {species_name}: {e}") if dataset.sizes["dim_electron_density"] <= 2 or dataset.sizes["dim_electron_temp"] <= 2: raise xr.AlignmentError(f"Alignment resulted in dataset sizes {dataset.sizes} for {species_name}.") + # Convert dimensionless coordinates back to physical quantities dataset["electron_density"] = dataset["dim_electron_density"] * reference_electron_density dataset["electron_temp"] = dataset["dim_electron_temp"] * reference_electron_temp + # 4. Standardize charge state indexing and attach global attributes dataset = align_rates_on_charge_states(dataset) - dataset = dataset.assign_attrs( atomic_number=config["species"][species_name]["atomic_number"], species_name=species_name, @@ -44,12 +54,10 @@ def read_rate_coeff(data_file_dir, species_name, config): created=datetime.date.today().strftime("%Y-%b-%d"), ) - dataset = write_global_attributes(dataset, config["globals"]) - - return dataset + return write_global_attributes(dataset, config["globals"]) def build_sorted_dictionary_of_rate_coefficients(config, species_name, data_file_dir): - """Parses configuration data to build a collection of rate coefficients sorted by year (most-recent first).""" + """Make a dictionary of rate coefficient datasets, ordered most-recent first.""" rate_coefficients = dict() years = dict() @@ -58,38 +66,31 @@ def build_sorted_dictionary_of_rate_coefficients(config, species_name, data_file config["data_file_config"], dataset_type ) + # Extract year for sorting; entries can be a simple int or [file, year] list if isinstance(file_to_read, int): years[dataset_type] = file_to_read elif isinstance(file_to_read, list) and len(file_to_read) == 2: years[dataset_type] = file_to_read[1] else: - raise NotImplementedError(f"Could not process entry: {file_to_read} for {species_name} {dataset_type}") + raise NotImplementedError(f"Unsupported config format for {species_name} {dataset_type}") if reader_key == "adf11": rate_dataset = build_adf11_rate_dataset( - data_file_dir, - species_name, - years[dataset_type], - dataset_type, - dataset_config, + data_file_dir, species_name, years[dataset_type], dataset_type, dataset_config, ) - else: - raise NotImplementedError( - f"No implementation for reading {reader_key} files." - ) + raise NotImplementedError(f"No implementation for reader: {reader_key}") rate_coefficients[dataset_type] = rate_dataset.rate_coefficient - # Sort the datasets so that the most recent data comes first - sorted_by_year = dict(sorted(years.items(), key=lambda item: item[1], reverse=True)) - - return {k: rate_coefficients[k] for k in sorted_by_year.keys()} + # Sort keys by year descending + sorted_keys = sorted(years, key=years.get, reverse=True) + return {k: rate_coefficients[k] for k in sorted_keys} def interpolate_rates_onto_matching_grids(config, rate_coefficients): - - # Since the rate coefficients are sorted by year, taking the first value - # gives the most recent data. + """Resample all rate coefficients to a uniform log-grid defined by the newest dataset.""" + + # Use the range of the most recent dataset to define the master grid most_recent_rate_coeff = list(rate_coefficients.values())[0] new_electron_density = np.logspace( @@ -105,16 +106,19 @@ def interpolate_rates_onto_matching_grids(config, rate_coefficients): ) interpolated_rate_coefficients = dict() - for key, value in rate_coefficients.items(): - interpolated_rate_coefficients[key] = value.groupby("dim_charge_state").map(interpolate_array, args=(new_electron_density, new_electron_temp)) + # Map interpolation across charge states + interpolated_rate_coefficients[key] = value.groupby("dim_charge_state").map( + interpolate_array, args=(new_electron_density, new_electron_temp) + ) return interpolated_rate_coefficients - def write_global_attributes(dataset: xr.Dataset, globals: dict) -> xr.Dataset: + """Attach global configuration parameters to the dataset as attributes or DataArrays.""" for attribute, value in globals.items(): if isinstance(value, dict): + # Complex attributes (with units) are added as coordinates/DataArrays if np.ndim(value["value"]) >= 1: dataset[attribute] = xr.DataArray( Quantity(value["value"], value["units"]), @@ -123,68 +127,58 @@ def write_global_attributes(dataset: xr.Dataset, globals: dict) -> xr.Dataset: else: dataset[attribute] = Quantity(value["value"], value["units"]) else: + # Simple metadata (strings/ints) added as attributes dataset[attribute] = value - return dataset - -def build_adf11_rate_dataset( - data_file_dir, species_name, year, dataset_type, dataset_config -): +def build_adf11_rate_dataset(data_file_dir, species_name, year, dataset_type, dataset_config): + """Read a specific ADF11 file and format it as a quantified xarray Dataset.""" from .adas_interface.read_adf11_file import read_adf11_file data = read_adf11_file(data_file_dir, species_name, year, dataset_type) - ds = xr.Dataset() - ds["species"] = species_name - ds["dataset"] = dataset_type - ds["charge"] = data["IZMAX"] - - electron_density = convert_units( - Quantity(10 ** data["DDENSD"][: data["IDMAXD"]], ureg.cm**-3), ureg.m**-3 - ) - electron_temp = Quantity(10 ** data["DTEVD"][: data["ITMAXD"]], ureg.eV) + # Log values stored in ADAS files are converted to linear scale if required + electron_density = convert_units(Quantity(10**data["DDENSD"][:data["IDMAXD"]], ureg.cm**-3), ureg.m**-3) + electron_temp = Quantity(10**data["DTEVD"][:data["ITMAXD"]], ureg.eV) - coefficient = data["DRCOFD"][: data["IZMAX"], : data["ITMAXD"], : data["IDMAXD"]] + coefficient = data["DRCOFD"][:data["IZMAX"], :data["ITMAXD"], :data["IDMAXD"]] if dataset_config["code"] <= 9: coefficient = 10**coefficient - input_units = dataset_config["stored_units"] - + # Create dimensionless indices for internal processing dim_electron_density = dimensionless_magnitude(electron_density / reference_electron_density) dim_electron_temp = dimensionless_magnitude(electron_temp / reference_electron_temp) - dim_charge_state = np.arange(data["IZMAX"]) - + rate_coefficient = xr.DataArray(coefficient, coords=dict( - dim_charge_state = dim_charge_state, + dim_charge_state = np.arange(data["IZMAX"]), dim_electron_temp = dim_electron_temp, dim_electron_density = dim_electron_density, - )).pint.quantify(input_units) + )).pint.quantify(dataset_config["stored_units"]) ds["rate_coefficient"] = convert_units(rate_coefficient, dataset_config["desired_units"]) - return ds - def align_rates_on_charge_states(dataset: xr.Dataset) -> xr.Dataset: - """For rates which are for k+1->k reactions, we shift these by one position - in the dim_charge_state dimension so that the kth position of a rate always - refers to the reactant species.""" - - dataset = dataset.pad( - pad_width=dict(dim_charge_state=(0, 1)), mode="constant", constant_values=0.0 - ) - dataset = dataset.assign_coords( - dim_charge_state=np.arange(dataset.sizes["dim_charge_state"]) - ) - - for key in [ + """ + Standardize charge state mapping so index 'k' always refers to the reactant. + + For k+1 -> k reactions (e.g. recombination), the rates are shifted so that + index k represents the species being recombined. + """ + # Pad to accommodate the N+1 charge state after shifting + dataset = dataset.pad(pad_width=dict(dim_charge_state=(0, 1)), mode="constant", constant_values=0.0) + dataset = dataset.assign_coords(dim_charge_state=np.arange(dataset.sizes["dim_charge_state"])) + + # Shift k+1 -> k processes + keys_to_shift = [ "effective_recombination", "charge_exchange_cross_coupling", "recombination_and_bremsstrahlung", "charge_exchange_emission", - ]: + ] + + for key in [k for k in keys_to_shift if k in dataset]: dataset[key] = dataset[key].roll(dim_charge_state=+1) - return dataset + return dataset \ No newline at end of file From 74cfa825539fc58e673681ace5e1782449d04dad Mon Sep 17 00:00:00 2001 From: Tom Body Date: Wed, 22 Apr 2026 14:27:39 -0400 Subject: [PATCH 19/25] Add new points for ne-tau --- radas/config.yaml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/radas/config.yaml b/radas/config.yaml index 4ccc42f..c777ff5 100644 --- a/radas/config.yaml +++ b/radas/config.yaml @@ -10,13 +10,13 @@ globals: # electron density (ne) * residence time (tau) (in m^-3 s) ne_tau: - value: [0.5e+17] + value: [0.5e+16, 0.5e+17, 0.5e+18] units: "m^-3 s" # Number of log-spaced points to use for electron density and electron temp. # Interpolation is used to map from raw data to points. - electron_density_resolution: 50 - electron_temp_resolution: 100 + electron_density_resolution: 20 + electron_temp_resolution: 80 data_file_config: adf11: From 9270a1bab6c1fddad3e0484aaa444d2390b6e554 Mon Sep 17 00:00:00 2001 From: Tom Body Date: Wed, 22 Apr 2026 14:37:10 -0400 Subject: [PATCH 20/25] Tidy up legend for plots --- radas/mavrin_reference/compare_to_mavrin.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/radas/mavrin_reference/compare_to_mavrin.py b/radas/mavrin_reference/compare_to_mavrin.py index f7c4974..7b8eda9 100644 --- a/radas/mavrin_reference/compare_to_mavrin.py +++ b/radas/mavrin_reference/compare_to_mavrin.py @@ -50,18 +50,21 @@ def compare_radas_to_mavrin_per_species(output_dir: Path, species: str, max_deca for i in range(ds.sizes["dim_ne_tau"]): ne_tau = ds.ne_tau.isel(dim_ne_tau=i).item() - Lz_radas.isel(dim_ne_tau=i).plot(ax=axs[0][0], label=f"{ne_tau:~P}", color=f"C{i}") + Lz_radas.isel(dim_ne_tau=i).plot(ax=axs[0][0], color=f"C{i}") if Lz_mavrin is not None: - Lz_mavrin.isel(dim_ne_tau=i).plot(ax=axs[0][1], color=f"C{i}") + Lz_mavrin.isel(dim_ne_tau=i).plot(ax=axs[0][1], color=f"C{i}", label=f"{ne_tau:~.1P}") + else: + axs[0][1].plot([], [], color=f"C{i}", label=f"{ne_tau:~.1P}") mean_charge_radas.isel(dim_ne_tau=i).plot(ax=axs[1][0], color=f"C{i}") if mean_charge_mavrin is not None: mean_charge_mavrin.isel(dim_ne_tau=i).plot(ax=axs[1][1], color=f"C{i}") - ds["coronal_Lz"].pint.to(ureg.W * ureg.m**3).plot(ax=axs[0][0], label="coronal", color="k") - ds["coronal_mean_charge_state"].pint.to(ureg.dimensionless).plot(ax=axs[1][0], label="coronal", color="k") + ds["coronal_Lz"].pint.to(ureg.W * ureg.m**3).plot(ax=axs[0][0], color="k") + axs[0][1].plot([], [], color="k", label="coronal") + ds["coronal_mean_charge_state"].pint.to(ureg.dimensionless).plot(ax=axs[1][0], color="k") - axs[0][0].legend() + axs[0][1].legend() Lz_radas_mag = magnitude_in_units(Lz_radas, ureg.W * ureg.m**3) Lz_coronal_mag = magnitude_in_units(ds["coronal_Lz"], ureg.W * ureg.m**3) @@ -94,4 +97,4 @@ def compare_radas_to_mavrin_per_species(output_dir: Path, species: str, max_deca if show: plt.show() - plt.savefig(output_dir / f"{species}.png") + plt.savefig(output_dir / f"{species}.png", dpi=300) From e1fd7679a88b5fe8558968fa2178d346de97695e Mon Sep 17 00:00:00 2001 From: Tom Body Date: Thu, 23 Apr 2026 12:14:16 -0400 Subject: [PATCH 21/25] Add reference electron density and temp back to datasets --- radas/read_rate_coeffs.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/radas/read_rate_coeffs.py b/radas/read_rate_coeffs.py index 82267bd..af618a0 100644 --- a/radas/read_rate_coeffs.py +++ b/radas/read_rate_coeffs.py @@ -43,6 +43,8 @@ def read_rate_coeff(data_file_dir, species_name, config): # Convert dimensionless coordinates back to physical quantities dataset["electron_density"] = dataset["dim_electron_density"] * reference_electron_density dataset["electron_temp"] = dataset["dim_electron_temp"] * reference_electron_temp + dataset["reference_electron_density"] = reference_electron_density + dataset["reference_electron_temp"] = reference_electron_temp # 4. Standardize charge state indexing and attach global attributes dataset = align_rates_on_charge_states(dataset) From 214f26d6eabd2b9c395f242f2a17ae031e1ce237 Mon Sep 17 00:00:00 2001 From: Tom Body Date: Thu, 23 Apr 2026 12:14:30 -0400 Subject: [PATCH 22/25] Remove git hash from datasets --- radas/read_rate_coeffs.py | 2 -- radas/shared.py | 13 ------------- 2 files changed, 15 deletions(-) diff --git a/radas/read_rate_coeffs.py b/radas/read_rate_coeffs.py index af618a0..15d7017 100644 --- a/radas/read_rate_coeffs.py +++ b/radas/read_rate_coeffs.py @@ -2,7 +2,6 @@ from .adas_interface.determine_adas_dataset_type import ( determine_reader_class_and_config, ) -from .shared import get_git_revision_short_hash from importlib.metadata import version, PackageNotFoundError import datetime import xarray as xr @@ -51,7 +50,6 @@ def read_rate_coeff(data_file_dir, species_name, config): dataset = dataset.assign_attrs( atomic_number=config["species"][species_name]["atomic_number"], species_name=species_name, - git_hash=get_git_revision_short_hash(), radas_version=radas_version, created=datetime.date.today().strftime("%Y-%b-%d"), ) diff --git a/radas/shared.py b/radas/shared.py index f8768b2..87ffa66 100644 --- a/radas/shared.py +++ b/radas/shared.py @@ -2,7 +2,6 @@ from pathlib import Path from importlib.resources import files -import subprocess import yaml default_config_file = files("radas").joinpath("config.yaml") @@ -11,18 +10,6 @@ library_extensions = [".a", ".so"] -def get_git_revision_short_hash() -> str: - try: - return ( - subprocess.check_output(["git", "rev-parse", "--short", "HEAD"]) - .decode("ascii") - .strip() - ) - except: # noqa:E722 - # If git isn't available (sometimes the case in tests), return a blank - return "UNDEFINED" - - def open_yaml_file(yaml_file: Path) -> dict: with open(yaml_file, "r") as file: return yaml.load(file, Loader=yaml.FullLoader) From 6454d912b4ed67121610b89720f38300b734f0fe Mon Sep 17 00:00:00 2001 From: Tom Body Date: Fri, 24 Apr 2026 18:24:02 -0400 Subject: [PATCH 23/25] Incorporate suggestions from @MishaVeldhoen --- .github/workflows/workflow_actions.yml | 6 ++++-- radas/cli.py | 27 +++++++++++++------------- 2 files changed, 18 insertions(+), 15 deletions(-) diff --git a/.github/workflows/workflow_actions.yml b/.github/workflows/workflow_actions.yml index 18abdc0..50252f1 100644 --- a/.github/workflows/workflow_actions.yml +++ b/.github/workflows/workflow_actions.yml @@ -10,6 +10,8 @@ on: env: FORCE_JAVASCRIPT_ACTIONS_TO_NODE24: true + # Use non-interactive matplotlib backend, to prevent the workflow from making + # interactive windows. MPLBACKEND: Agg concurrency: @@ -57,7 +59,7 @@ jobs: test_radas: name: Test radas command runs-on: ubuntu-latest - if: (startsWith(github.ref, 'refs/tags/') || github.event_name == 'workflow_dispatch') + if: (github.event_name == 'release' || github.event_name == 'workflow_dispatch') permissions: contents: write @@ -81,7 +83,7 @@ jobs: - name: Create and Upload Release Asset # Only upload if this is a tag (skip for manual workflow_dispatch runs) - if: startsWith(github.ref, 'refs/tags/') + if: github.event_name == 'release' env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} run: | diff --git a/radas/cli.py b/radas/cli.py index f501043..83e8146 100644 --- a/radas/cli.py +++ b/radas/cli.py @@ -70,21 +70,22 @@ def run_radas_cli( debug=debug, ) - # nullcontext does nothing. - debug_context = contextlib.nullcontext() - if debug: - try: - # If --debug and ipdb is installed, switch - # to use the launch_ipdb_on_exception context - from ipdb import launch_ipdb_on_exception - debug_context = launch_ipdb_on_exception() - except ModuleNotFoundError: - pass - - with debug_context: + with _post_mortem_debugger(): + run_radas(**kwargs) + else: run_radas(**kwargs) +def _post_mortem_debugger(): + """Context manager that drops into ipdb on unhandled exceptions, or a no-op if ipdb is not installed.""" + try: + from ipdb import launch_ipdb_on_exception + except ModuleNotFoundError: + print("Warning: --debug set but ipdb is not installed; " + "install with `pip install ipdb` for post-mortem debugging.") + return contextlib.nullcontext() + return launch_ipdb_on_exception() + def run_radas( directory: Path, config: Optional[str], @@ -133,7 +134,7 @@ def run_radas( (species_name in species) or (species == ("all",)) ): datasets[species_name] = read_rate_coeff( - data_file_dir, species_name, configuration + data_file_dir, species_name, configuration, debug=debug ) output_dir.mkdir(exist_ok=True, parents=True) From b295b3f916ff8fdc58ea8269717e8977f0e0c5ab Mon Sep 17 00:00:00 2001 From: Tom Body Date: Fri, 24 Apr 2026 18:27:32 -0400 Subject: [PATCH 24/25] Add a warning for off-grid extrapolation (switch on with --debug) --- radas/interpolate_rates.py | 35 +++++++++++++++++++++++++++++++++++ radas/read_rate_coeffs.py | 27 ++++++++++++++++++++------- 2 files changed, 55 insertions(+), 7 deletions(-) diff --git a/radas/interpolate_rates.py b/radas/interpolate_rates.py index 2af7c34..88629cb 100644 --- a/radas/interpolate_rates.py +++ b/radas/interpolate_rates.py @@ -3,6 +3,13 @@ import numpy as np from scipy.interpolate import RectBivariateSpline from numpy.typing import NDArray +import warnings + +def is_significantly_below(requested, limit): + return requested < limit and not np.isclose(requested, limit) + +def is_significantly_above(requested, limit): + return requested > limit and not np.isclose(requested, limit) def interpolate_array( array: xr.DataArray, @@ -28,6 +35,34 @@ def interpolate_array( if np.any(array <= 0.0): raise NotImplementedError("Cannot log-interpolate rate coefficients containing zeros.") + # Check if extrapolation is needed and raise a warning if this is the case. + out_of_bounds_msg = [] + + req_dens_min, req_dens_max = new_electron_density.min(), new_electron_density.max() + grid_dens_min, grid_dens_max = array.dim_electron_density.min(), array.dim_electron_density.max() + + if is_significantly_below(req_dens_min, grid_dens_min) or is_significantly_above(req_dens_max, grid_dens_max): + out_of_bounds_msg.append( + f"Density requested [{req_dens_min:.2e}, {req_dens_max:.2e}] " + f"exceeds grid [{grid_dens_min:.2e}, {grid_dens_max:.2e}]." + ) + + # Check Temperature Bounds + req_temp_min, req_temp_max = new_electron_temp.min(), new_electron_temp.max() + grid_temp_min, grid_temp_max = array.dim_electron_temp.min(), array.dim_electron_temp.max() + + if is_significantly_below(req_temp_min, grid_temp_min) or is_significantly_above(req_temp_max, grid_temp_max): + out_of_bounds_msg.append( + f"Temperature requested [{req_temp_min:.2e}, {req_temp_max:.2e}] " + f"exceeds grid [{grid_temp_min:.2e}, {grid_temp_max:.2e}]." + ) + + if out_of_bounds_msg: + full_msg = "Nearest-neighbour extrapolation used for off-grid values: " + " ".join(out_of_bounds_msg) + warnings.warn(full_msg, RuntimeWarning) + + # ------------------------------------ + # Prepare original grid and data in log10 space x = np.log10(array.dim_electron_density) y = np.log10(array.dim_electron_temp) diff --git a/radas/read_rate_coeffs.py b/radas/read_rate_coeffs.py index 15d7017..c0dde01 100644 --- a/radas/read_rate_coeffs.py +++ b/radas/read_rate_coeffs.py @@ -6,13 +6,14 @@ import datetime import xarray as xr import numpy as np +import warnings from .interpolate_rates import interpolate_array # Reference units for non-dimensionalizing coordinates reference_electron_density = Quantity(1.0, ureg.m**-3) reference_electron_temp = Quantity(1.0, ureg.eV) -def read_rate_coeff(data_file_dir, species_name, config): +def read_rate_coeff(data_file_dir, species_name, config, debug=False): """ Main pipeline to assemble an atomic rate dataset for a specific species. @@ -28,7 +29,7 @@ def read_rate_coeff(data_file_dir, species_name, config): rate_coefficients = build_sorted_dictionary_of_rate_coefficients(config, species_name, data_file_dir) # 2. Resample all datasets to a common resolution - rate_coefficients = interpolate_rates_onto_matching_grids(config, rate_coefficients) + rate_coefficients = interpolate_rates_onto_matching_grids(config, species_name, rate_coefficients, debug) # 3. Merge individual datasets (e.g., recombination, ionization) into one try: @@ -87,7 +88,7 @@ def build_sorted_dictionary_of_rate_coefficients(config, species_name, data_file sorted_keys = sorted(years, key=years.get, reverse=True) return {k: rate_coefficients[k] for k in sorted_keys} -def interpolate_rates_onto_matching_grids(config, rate_coefficients): +def interpolate_rates_onto_matching_grids(config, species_name, rate_coefficients, debug): """Resample all rate coefficients to a uniform log-grid defined by the newest dataset.""" # Use the range of the most recent dataset to define the master grid @@ -107,10 +108,22 @@ def interpolate_rates_onto_matching_grids(config, rate_coefficients): interpolated_rate_coefficients = dict() for key, value in rate_coefficients.items(): - # Map interpolation across charge states - interpolated_rate_coefficients[key] = value.groupby("dim_charge_state").map( - interpolate_array, args=(new_electron_density, new_electron_temp) - ) + with warnings.catch_warnings(record=True) as captured_warnings: + warnings.simplefilter("always") + + # Map interpolation across charge states + interpolated_rate_coefficients[key] = value.groupby("dim_charge_state").map( + interpolate_array, args=(new_electron_density, new_electron_temp) + ) + + if debug: + for w in captured_warnings: + warnings.warn_explicit( + message=f"when interpolating {key} for {species_name}. {w.message}", + category=w.category, + filename=w.filename, + lineno=w.lineno, + ) return interpolated_rate_coefficients From 844560d9b193dba5b33091da6d4a1af63b243535 Mon Sep 17 00:00:00 2001 From: Tom Body Date: Mon, 27 Apr 2026 13:16:53 -0400 Subject: [PATCH 25/25] Incorporate comments from @MishaVeldhoen --- .github/workflows/workflow_actions.yml | 25 +++++-------------------- radas/cli.py | 2 +- radas/read_rate_coeffs.py | 8 ++++---- 3 files changed, 10 insertions(+), 25 deletions(-) diff --git a/.github/workflows/workflow_actions.yml b/.github/workflows/workflow_actions.yml index 50252f1..d3f439f 100644 --- a/.github/workflows/workflow_actions.yml +++ b/.github/workflows/workflow_actions.yml @@ -10,8 +10,8 @@ on: env: FORCE_JAVASCRIPT_ACTIONS_TO_NODE24: true - # Use non-interactive matplotlib backend, to prevent the workflow from making - # interactive windows. + # Use non-interactive matplotlib backend, to prevent the workflow from + # trying to make interactive windows. MPLBACKEND: Agg concurrency: @@ -55,7 +55,7 @@ jobs: radas_config -o ./new_config.yaml radas -s hydrogen -c ./new_config.yaml - # --- JOB 2: RUN RADAS COMMAND --- + # --- JOB 2: RUN RADAS COMMAND, CHECK IT WORKS WITH DEFAULT CONFIG --- test_radas: name: Test radas command runs-on: ubuntu-latest @@ -81,26 +81,11 @@ jobs: - name: Run radas command run: poetry run radas -vvv - - name: Create and Upload Release Asset - # Only upload if this is a tag (skip for manual workflow_dispatch runs) - if: github.event_name == 'release' - env: - GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} - run: | - if [ -d "./radas_dir" ]; then - zip -r radas_dir.zip ./radas_dir - gh release create ${{ github.ref_name }} --generate-notes || true - gh release upload ${{ github.ref_name }} radas_dir.zip --clobber - else - echo "::error::radas_dir not found! Skipping upload." - exit 1 - fi - # --- JOB 3: BUILD RELEASE ARTIFACTS --- build_release: name: Build Release needs: [pytest, test_radas] - if: startsWith(github.ref, 'refs/tags') + if: github.event_name == 'release' runs-on: ubuntu-latest steps: - uses: actions/checkout@v5 @@ -122,7 +107,7 @@ jobs: publish: name: Publish to PyPI needs: build_release - if: startsWith(github.ref, 'refs/tags') + if: github.event_name == 'release' runs-on: ubuntu-latest environment: name: pypi-publish diff --git a/radas/cli.py b/radas/cli.py index 83e8146..e05f349 100644 --- a/radas/cli.py +++ b/radas/cli.py @@ -134,7 +134,7 @@ def run_radas( (species_name in species) or (species == ("all",)) ): datasets[species_name] = read_rate_coeff( - data_file_dir, species_name, configuration, debug=debug + data_file_dir, species_name, configuration, verbose=verbose, ) output_dir.mkdir(exist_ok=True, parents=True) diff --git a/radas/read_rate_coeffs.py b/radas/read_rate_coeffs.py index c0dde01..a3e279a 100644 --- a/radas/read_rate_coeffs.py +++ b/radas/read_rate_coeffs.py @@ -13,7 +13,7 @@ reference_electron_density = Quantity(1.0, ureg.m**-3) reference_electron_temp = Quantity(1.0, ureg.eV) -def read_rate_coeff(data_file_dir, species_name, config, debug=False): +def read_rate_coeff(data_file_dir, species_name, config, verbose=0): """ Main pipeline to assemble an atomic rate dataset for a specific species. @@ -29,7 +29,7 @@ def read_rate_coeff(data_file_dir, species_name, config, debug=False): rate_coefficients = build_sorted_dictionary_of_rate_coefficients(config, species_name, data_file_dir) # 2. Resample all datasets to a common resolution - rate_coefficients = interpolate_rates_onto_matching_grids(config, species_name, rate_coefficients, debug) + rate_coefficients = interpolate_rates_onto_matching_grids(config, species_name, rate_coefficients, verbose=verbose) # 3. Merge individual datasets (e.g., recombination, ionization) into one try: @@ -88,7 +88,7 @@ def build_sorted_dictionary_of_rate_coefficients(config, species_name, data_file sorted_keys = sorted(years, key=years.get, reverse=True) return {k: rate_coefficients[k] for k in sorted_keys} -def interpolate_rates_onto_matching_grids(config, species_name, rate_coefficients, debug): +def interpolate_rates_onto_matching_grids(config, species_name, rate_coefficients, verbose): """Resample all rate coefficients to a uniform log-grid defined by the newest dataset.""" # Use the range of the most recent dataset to define the master grid @@ -116,7 +116,7 @@ def interpolate_rates_onto_matching_grids(config, species_name, rate_coefficient interpolate_array, args=(new_electron_density, new_electron_temp) ) - if debug: + if verbose: for w in captured_warnings: warnings.warn_explicit( message=f"when interpolating {key} for {species_name}. {w.message}",