From c1c1ff50ce5d1e3a98df16bd988eca02988ce4bf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jakub=20Kali=C5=84ski?= Date: Wed, 15 Jul 2026 16:49:22 +0200 Subject: [PATCH] Add prototype of easier comparison with experimental data --- requirements.txt | 1 + validation/README.md | 303 +++++++++++++++++++++++++++++++++ validation/__init__.py | 15 ++ validation/comparison.py | 254 +++++++++++++++++++++++++++ validation/config.py | 142 +++++++++++++++ validation/config.yaml | 23 +++ validation/config_example.yaml | 57 +++++++ validation/data_utils.py | 209 +++++++++++++++++++++++ validation/plots.py | 236 +++++++++++++++++++++++++ validation/run_comparison.py | 90 ++++++++++ 10 files changed, 1330 insertions(+) create mode 100644 validation/README.md create mode 100644 validation/__init__.py create mode 100644 validation/comparison.py create mode 100644 validation/config.py create mode 100644 validation/config.yaml create mode 100644 validation/config_example.yaml create mode 100644 validation/data_utils.py create mode 100644 validation/plots.py create mode 100644 validation/run_comparison.py diff --git a/requirements.txt b/requirements.txt index 4f7f27b..5197d54 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,6 +6,7 @@ mpmath pandas seaborn numba +pyyaml isort flake8 flake8-isort \ No newline at end of file diff --git a/validation/README.md b/validation/README.md new file mode 100644 index 0000000..112d21e --- /dev/null +++ b/validation/README.md @@ -0,0 +1,303 @@ +# Experimental Data Comparison Module + +This module provides tools for comparing IonTracks simulation results with experimental measurements of ion recombination in ionization chambers. + +## Overview + +The validation module allows you to: +- Load experimental data from CSV files +- Run IonTracks simulations with configurable parameters (backend, grid size, etc.) +- Compare experimental and simulated results for both initial and general recombination +- Generate publication-quality plots and analysis reports + +## Installation + +The validation module is part of the IonTracks-Cython package. All required dependencies are included in the project's `requirements.txt`. Install them with: + +```bash +pip install -r requirements.txt +``` + +This will install all necessary packages including pandas, numpy, matplotlib, seaborn, and pyyaml. + +## Quick Start + +### 1. Create a Configuration File + +First, create a configuration template. It's recommended to place it in the same directory as your experimental data or in the `validation/` directory: + +```bash +# Option 1: In validation directory (good for project-wide configs) +python -m validation.run_comparison --create-template validation/config.yaml + +# Option 2: In experimental data directory (good for experiment-specific configs) +python -m validation.run_comparison --create-template experimental-data/2019_DCPT/config.yaml +``` + +This creates a template YAML file that you can edit to specify your experimental data path and simulation parameters. + +### 2. Edit Configuration + +Edit `config.yaml` to specify: + +```yaml +experimental_data_path: "path/to/your/experimental_data.csv" +output_dir: "results" + +simulation: + backend: "cython" # Options: cython, python, numba, cupy, parallel + voltage_V: 200.0 + electrode_gap_cm: 0.2 + particle: "proton" + RDD_model: "Gauss" # Options: Gauss, Geiss + grid_size_um: 5.0 + a0_nm: 8.0 + use_beta: False + seed: null # null for random seed + +compare_initial_recombination: true +compare_continuous_beam: true +``` + +### 3. Prepare Experimental Data + +Your experimental data CSV file should contain the following columns (exact names are flexible, see below): + +**Required columns:** +- **Energy** (MeV): Proton or ion energy + - Acceptable names: `Energy_MeV`, `energy_MeV`, `E_MeV`, `energy` + +- **k_s**: Recombination correction factor + - Acceptable names: `k_s`, `ks`, `recombination_factor`, `collection_efficiency` + +- **Dose rate** (Gy/s): Dose rate in air or water + - For air: `dose_rate_air_Gy_s`, `dose_rate_air`, `doserate_air_Gy_s` + - For water: `dose_rate_water_Gy_s`, `dose_rate_water`, `doserate_water_Gy_s` + - **Important**: Use air dose rate if your ionization chamber is filled with air (recommended) + +**Optional columns:** +- Any other metadata columns will be preserved in output files + +**Example CSV structure:** + +```csv +Energy_MeV,dose_rate_air_Gy_s,k_s +150,3.665,1.0025 +150,1.866,1.0013 +226,5.61,1.0038 +226,2.877,1.0019 +``` + +**Note on dose rate units:** +- The code automatically converts Gy/s to Gy/min internally for calculations +- Always use the dose rate that matches your ionization chamber medium (air vs. water) +- Using the wrong medium (e.g., water dose rate for air-filled chamber) introduces ~13-14% error + +### 4. Run Comparison + +```bash +# Use the path to your config file +python -m validation.run_comparison validation/config.yaml +# or +python -m validation.run_comparison experimental-data/2019_DCPT/config.yaml +``` + +This will: +1. Load and validate your experimental data +2. Run IonTracks simulations for initial recombination (if enabled) +3. Run IonTracks simulations for continuous beam/general recombination (if enabled) +4. Compare results and calculate errors +5. Generate plots +6. Save results to CSV files + +## Configuration Options + +### Simulation Parameters + +#### Backend Selection + +Choose the computational backend: + +- **`cython`** (default): Fast, compiled Cython code +- **`python`**: Pure Python (slower, but no compilation needed) +- **`numba`**: JIT-compiled NumPy code +- **`cupy`**: GPU-accelerated (requires CUDA-capable GPU) +- **`parallel`**: Parallel NumPy implementation + +#### Physics Parameters + +- **`voltage_V`**: Applied voltage in volts (typically 200-300 V) +- **`electrode_gap_cm`**: Electrode gap in cm (typically 0.1-0.2 cm) +- **`particle`**: Particle type (e.g., `"proton"`, `"carbon"`, `"helium"`) +- **`RDD_model`**: Radial dose distribution model + - `"Gauss"`: Gaussian distribution (default) + - `"Geiss"`: Geiss model +- **`grid_size_um`**: Spatial grid resolution in micrometers (smaller = more accurate but slower) +- **`a0_nm`**: Track core radius parameter in nanometers +- **`use_beta`**: Scale track core by β = v/c (usually False) + +#### Computational Options + +- **`seed`**: Random seed for reproducibility (null for random) +- **`SHOW_PLOT`**: Show simulation plots during calculation (False) +- **`PRINT_parameters`**: Print detailed parameters (False) +- **`debug`**: Enable debug mode (False) + +### Comparison Options + +- **`compare_initial_recombination`**: Compare initial recombination (fast, recommended) +- **`compare_continuous_beam`**: Compare general recombination (slow, may take minutes to hours) +- **`dose_rate_column`**: Which dose rate column to use (auto-detected if not specified) +- **`energy_column`**: Which energy column to use (auto-detected if not specified) +- **`ks_column`**: Which k_s column to use (auto-detected if not specified) + +### Plotting Options + +- **`plot_format`**: Output format (`"png"`, `"pdf"`, `"svg"`) +- **`plot_dpi`**: Resolution for raster formats (300 recommended) +- **`plot_style`**: Matplotlib style (`"whitegrid"`, `"darkgrid"`, etc.) + +## Output Files + +The comparison generates the following files in the output directory: + +### CSV Files + +1. **`comparison_initial.csv`**: Initial recombination comparison + - Columns: `Energy_MeV`, `k_s_experimental`, `k_s_IonTracks`, `difference`, `relative_error_%` + +2. **`comparison_continuous.csv`**: Full continuous beam comparison + - Includes all original experimental columns plus IonTracks results and errors + +3. **`comparison_continuous_clean.csv`**: Simplified continuous beam comparison + - Only essential columns: energy, dose rate, k_s values, and errors + +### Plot Files + +1. **`initial_recombination_comparison.png`**: Initial recombination vs. energy +2. **`continuous_beam_comparison.png`**: Continuous beam results vs. dose rate (all energies) +3. **`relative_error.png`**: Relative error as function of dose rate + +## Understanding Results + +### Initial Recombination + +Initial recombination occurs at very low dose rates where individual ion tracks don't overlap. The comparison: +- Extracts the minimum k_s value for each energy (lowest dose rate = initial recombination) +- Compares with IonTracks single-track calculations +- Shows how well the model predicts initial recombination + +### General Recombination (Continuous Beam) + +General recombination occurs at higher dose rates where track overlap becomes significant. The comparison: +- Runs IonTracks for each experimental condition (energy + dose rate) +- Compares simulated and experimental k_s values +- Calculates relative errors to assess model accuracy + +### Error Metrics + +- **Absolute error**: `k_s_IonTracks - k_s_experimental` +- **Relative error**: `100 × (k_s_IonTracks - k_s_experimental) / k_s_experimental` +- Statistics include mean, median, min, max, and standard deviation + +## Example Workflow + +### Example 1: Quick Initial Recombination Check + +```yaml +experimental_data_path: "experimental-data/2019_DCPT/recombination_200V_data.csv" +output_dir: "results/initial_only" +simulation: + backend: "cython" + voltage_V: 200.0 + electrode_gap_cm: 0.2 +compare_initial_recombination: true +compare_continuous_beam: false +``` + +### Example 2: Full Comparison with Custom Parameters + +```yaml +experimental_data_path: "my_data.csv" +output_dir: "results/custom" +simulation: + backend: "numba" + voltage_V: 300.0 + electrode_gap_cm: 0.1 + particle: "proton" + RDD_model: "Geiss" + grid_size_um: 3.0 + seed: 42 +compare_initial_recombination: true +compare_continuous_beam: true +``` + +### Example 3: High-Resolution Simulation + +```yaml +simulation: + backend: "cython" + grid_size_um: 2.0 # Smaller = more accurate but slower + a0_nm: 10.0 +``` + +## Troubleshooting + +### Common Issues + +1. **"Could not find energy column"** + - Check that your CSV has a column named one of: `Energy_MeV`, `energy_MeV`, `E_MeV`, `energy` + - Or specify the exact column name in config: `energy_column: "your_column_name"` + +2. **"Failed to calculate any results"** + - Check that parameters are physically reasonable + - Try a different backend (e.g., `python` instead of `cython`) + - Enable debug mode: `simulation.debug: true` + +3. **"Slow performance"** + - Use `cython` or `numba` backend instead of `python` + - Increase `grid_size_um` (less accurate but faster) + - Disable `compare_continuous_beam` for quick initial checks + +4. **"Import errors"** + - Ensure IonTracks-Cython package is properly installed + - Check that you're running from the correct directory + - Verify all dependencies are installed + +## Programmatic Usage + +You can also use the validation module programmatically: + +```python +from pathlib import Path +from validation.config import load_config +from validation.comparison import ComparisonRunner +from validation.plots import generate_all_plots + +# Load configuration +config = load_config("config.yaml") + +# Run comparison +runner = ComparisonRunner(config) +runner.run_all() +runner.save_results() + +# Generate plots +generate_all_plots( + runner.comparison_initial, + runner.comparison_continuous, + config, +) +``` + +## References + +For more information about IonTracks and validation studies: + +- Christensen, J.B. et al. (2020) "Mapping initial and general recombination in scanning proton pencil beams" *Phys. Med. Biol.* 65 115003 +- Christensen, J.B., Tölli, H., Bassler, N. (2016) "A general algorithm for calculation of recombination losses in ionization chambers exposed to ion beams" *Medical Physics* 43.10: 5484-92 + +## License + +Part of the IonTracks-Cython project. + diff --git a/validation/__init__.py b/validation/__init__.py new file mode 100644 index 0000000..75c6610 --- /dev/null +++ b/validation/__init__.py @@ -0,0 +1,15 @@ +""" +Validation module for comparing IonTracks simulation results with experimental data. + +This module provides tools for: +- Loading and validating experimental data +- Running IonTracks simulations with configurable parameters +- Comparing experimental and simulated results +- Generating plots and analysis reports +""" + +from validation.comparison import ComparisonRunner +from validation.data_utils import load_experimental_data, validate_experimental_data + +__all__ = ["ComparisonRunner", "load_experimental_data", "validate_experimental_data"] + diff --git a/validation/comparison.py b/validation/comparison.py new file mode 100644 index 0000000..1a5771d --- /dev/null +++ b/validation/comparison.py @@ -0,0 +1,254 @@ +""" +Main comparison logic for experimental data vs. IonTracks simulations. +""" + +import sys +from pathlib import Path +from typing import Optional + +import numpy as np +import pandas as pd + +# Add parent directory to path for imports +sys.path.insert(0, str(Path(__file__).parent.parent)) + +from hadrons.functions import ks_initial_IonTracks, IonTracks_continuous_beam + +from validation.config import ComparisonConfig, SimulationConfig +from validation.data_utils import load_experimental_data + + +class ComparisonRunner: + """Runner for comparing experimental data with IonTracks simulations.""" + + def __init__(self, config: ComparisonConfig): + """ + Initialize comparison runner. + + Parameters + ---------- + config : ComparisonConfig + Configuration object + """ + self.config = config + self.config.output_dir.mkdir(parents=True, exist_ok=True) + + # Load experimental data + self.exp_data = load_experimental_data( + config.experimental_data_path, + energy_column=config.energy_column, + ks_column=config.ks_column, + dose_rate_column=config.dose_rate_column, + ) + + # Results storage + self.initial_results: Optional[pd.DataFrame] = None + self.continuous_results: Optional[pd.DataFrame] = None + self.comparison_initial: Optional[pd.DataFrame] = None + self.comparison_continuous: Optional[pd.DataFrame] = None + + def run_initial_recombination(self) -> pd.DataFrame: + """ + Run initial recombination calculations and compare with experimental data. + + Returns + ------- + pd.DataFrame + Comparison results DataFrame + """ + print("\n" + "=" * 70) + print("INITIAL RECOMBINATION COMPARISON") + print("=" * 70) + + sim_config = self.config.simulation + energies_MeV = sorted(self.exp_data["Energy_MeV"].unique()) + + initial_results_list = [] + + for energy_MeV in energies_MeV: + print(f"Calculating for {energy_MeV} MeV...") + try: + result = ks_initial_IonTracks( + E_MeV_u=energy_MeV, + voltage_V=sim_config.voltage_V, + electrode_gap_cm=sim_config.electrode_gap_cm, + particle=sim_config.particle, + RDD_model=sim_config.RDD_model, + grid_size_um=sim_config.grid_size_um, + a0_nm=sim_config.a0_nm, + use_beta=sim_config.use_beta, + debug=sim_config.debug, + SHOW_PLOT=sim_config.SHOW_PLOT, + ) + initial_results_list.append(result) + except Exception as e: + print(f" Error for {energy_MeV} MeV: {e}") + continue + + if not initial_results_list: + raise RuntimeError("Failed to calculate any results for Initial Recombination") + + self.initial_results = pd.concat(initial_results_list, ignore_index=True) + + # Compare with experimental data (lowest dose rate = initial recombination) + comparison_list = [] + for energy_MeV in energies_MeV: + exp_subset = self.exp_data[self.exp_data["Energy_MeV"] == energy_MeV] + k_s_exp_min = exp_subset["k_s"].min() # Lowest dose rate + + iontracks_subset = self.initial_results[ + self.initial_results["E_MeV_u"] == energy_MeV + ] + if len(iontracks_subset) > 0: + k_s_iontracks = iontracks_subset["ks"].values[0] + comparison_list.append( + { + "Energy_MeV": energy_MeV, + "k_s_experimental": k_s_exp_min, + "k_s_IonTracks": k_s_iontracks, + "difference": k_s_iontracks - k_s_exp_min, + "relative_error_%": 100 + * (k_s_iontracks - k_s_exp_min) + / k_s_exp_min, + } + ) + + self.comparison_initial = pd.DataFrame(comparison_list) + + print("\nInitial Recombination Comparison:") + print(self.comparison_initial.to_string(index=False)) + + return self.comparison_initial + + def run_continuous_beam(self) -> pd.DataFrame: + """ + Run continuous beam calculations and compare with experimental data. + + Returns + ------- + pd.DataFrame + Comparison results DataFrame + """ + print("\n" + "=" * 70) + print("CONTINUOUS BEAM COMPARISON (General Recombination)") + print("=" * 70) + print("WARNING: This may take several minutes...") + + sim_config = self.config.simulation + continuous_results_list = [] + + # Determine seed + if sim_config.seed is None: + import random + seed = random.randint(1, int(1e7)) + else: + seed = sim_config.seed + + total_rows = len(self.exp_data) + for idx, row in self.exp_data.iterrows(): + if idx % max(1, total_rows // 10) == 0: + progress = 100 * (idx + 1) / total_rows + print(f"Progress: {idx+1}/{total_rows} ({progress:.1f}%)") + + try: + result = IonTracks_continuous_beam( + E_MeV_u=row["Energy_MeV"], + voltage_V=sim_config.voltage_V, + doserate_Gy_min=row["dose_rate_Gy_min"], + electrode_gap_cm=sim_config.electrode_gap_cm, + particle=sim_config.particle, + grid_size_um=sim_config.grid_size_um, + backend=sim_config.backend, + PRINT_parameters=sim_config.PRINT_parameters, + SHOW_PLOT=sim_config.SHOW_PLOT, + myseed=seed + idx, # Different seed for each calculation + ) + continuous_results_list.append(result) + except Exception as e: + print(f" Error for row {idx}: {e}") + # Add NaN as result to preserve indices + continuous_results_list.append( + pd.DataFrame([{"ks_IonTracks": np.nan, "E_MeV_u": row["Energy_MeV"]}]) + ) + + if not continuous_results_list: + raise RuntimeError("Failed to calculate any results for Continuous Beam") + + continuous_df = pd.concat(continuous_results_list, ignore_index=True) + self.continuous_results = continuous_df + + # Create comparison DataFrame + comparison_continuous = self.exp_data.copy() + comparison_continuous["k_s_IonTracks"] = continuous_df["ks_IonTracks"].values + comparison_continuous["difference"] = ( + comparison_continuous["k_s_IonTracks"] - comparison_continuous["k_s"] + ) + comparison_continuous["relative_error_%"] = ( + 100 + * comparison_continuous["difference"] + / comparison_continuous["k_s"] + ) + comparison_continuous["absolute_error"] = abs( + comparison_continuous["difference"] + ) + + # Remove rows with NaN + comparison_continuous = comparison_continuous.dropna(subset=["k_s_IonTracks"]) + self.comparison_continuous = comparison_continuous + + # Print summary + print("\nContinuous Beam Comparison Summary:") + print(f"Number of successful calculations: {len(comparison_continuous)}/{len(self.exp_data)}") + if len(comparison_continuous) > 0: + print(f"Mean absolute error: {comparison_continuous['absolute_error'].mean():.6f}") + print(f"Mean relative error: {comparison_continuous['relative_error_%'].mean():.3f}%") + print(f"Median relative error: {comparison_continuous['relative_error_%'].median():.3f}%") + print(f"Max relative error: {comparison_continuous['relative_error_%'].max():.3f}%") + print(f"Min relative error: {comparison_continuous['relative_error_%'].min():.3f}%") + print(f"Standard deviation: {comparison_continuous['relative_error_%'].std():.3f}%") + + return self.comparison_continuous + + def run_all(self): + """Run both initial recombination and continuous beam comparisons.""" + if self.config.compare_initial_recombination: + self.run_initial_recombination() + + if self.config.compare_continuous_beam: + self.run_continuous_beam() + + def save_results(self): + """Save comparison results to CSV files.""" + print("\n" + "=" * 70) + print("SAVING RESULTS") + print("=" * 70) + + if self.comparison_initial is not None and len(self.comparison_initial) > 0: + output_path = self.config.output_dir / "comparison_initial.csv" + self.comparison_initial.to_csv(output_path, index=False) + print(f"✓ Saved: {output_path}") + + if self.comparison_continuous is not None and len(self.comparison_continuous) > 0: + # Save full comparison + output_path = self.config.output_dir / "comparison_continuous.csv" + self.comparison_continuous.to_csv(output_path, index=False) + print(f"✓ Saved: {output_path}") + + # Save clean version + clean_cols = [ + "Energy_MeV", + "dose_rate_Gy_s", + "k_s", + "k_s_IonTracks", + "difference", + "relative_error_%", + "absolute_error", + ] + available_cols = [col for col in clean_cols if col in self.comparison_continuous.columns] + comparison_clean = self.comparison_continuous[available_cols] + output_path_clean = self.config.output_dir / "comparison_continuous_clean.csv" + comparison_clean.to_csv(output_path_clean, index=False) + print(f"✓ Saved: {output_path_clean}") + + print(f"\nAll results saved in: {self.config.output_dir}") + diff --git a/validation/config.py b/validation/config.py new file mode 100644 index 0000000..0c2540d --- /dev/null +++ b/validation/config.py @@ -0,0 +1,142 @@ +""" +Configuration management for experimental data comparison. + +Supports loading configuration from YAML files with validation and defaults. +""" + +from dataclasses import dataclass, field +from pathlib import Path +from typing import Literal, Optional + +import yaml + + +@dataclass +class SimulationConfig: + """Configuration for IonTracks simulations.""" + + # Backend selection + backend: Literal["cython", "python", "numba", "cupy", "parallel"] = "cython" + + # Experimental setup parameters + voltage_V: float = 200.0 + electrode_gap_cm: float = 0.2 + particle: str = "proton" + + # Simulation parameters + RDD_model: Literal["Gauss", "Geiss"] = "Gauss" + grid_size_um: float = 5.0 + a0_nm: float = 8.0 + use_beta: bool = False + seed: Optional[int] = None + + # Optional simulation flags + SHOW_PLOT: bool = False + PRINT_parameters: bool = False + debug: bool = False + + +@dataclass +class ComparisonConfig: + """Configuration for comparison analysis.""" + + # Paths + experimental_data_path: Path + output_dir: Path = field(default_factory=lambda: Path("results")) + + # Simulation configuration + simulation: SimulationConfig = field(default_factory=SimulationConfig) + + # Comparison options + compare_initial_recombination: bool = True + compare_continuous_beam: bool = True + + # Data processing options (None for auto-detection) + dose_rate_column: Optional[str] = None # Auto-detect if None + energy_column: Optional[str] = None # Auto-detect if None + ks_column: Optional[str] = None # Auto-detect if None + + # Plotting options + plot_format: str = "png" + plot_dpi: int = 300 + plot_style: str = "whitegrid" + + def __post_init__(self): + """Convert string paths to Path objects if needed.""" + if isinstance(self.experimental_data_path, str): + self.experimental_data_path = Path(self.experimental_data_path) + if isinstance(self.output_dir, str): + self.output_dir = Path(self.output_dir) + + +def load_config(config_path: Path) -> ComparisonConfig: + """ + Load configuration from YAML file. + + Parameters + ---------- + config_path : Path + Path to YAML configuration file + + Returns + ------- + ComparisonConfig + Loaded configuration object + """ + config_path = Path(config_path) + + if not config_path.exists(): + raise FileNotFoundError(f"Configuration file not found: {config_path}") + + with open(config_path, "r") as f: + config_dict = yaml.safe_load(f) + + # Extract simulation config + sim_dict = config_dict.get("simulation", {}) + simulation_config = SimulationConfig(**sim_dict) + + # Create comparison config + config_dict["simulation"] = simulation_config + config = ComparisonConfig(**config_dict) + + return config + + +def save_config_template(template_path: Path): + """ + Save a template configuration file. + + Parameters + ---------- + template_path : Path + Path where to save the template + """ + template = { + "experimental_data_path": "path/to/experimental_data.csv", + "output_dir": "results", + "simulation": { + "backend": "cython", # Options: cython, python, numba, cupy, parallel + "voltage_V": 200.0, + "electrode_gap_cm": 0.2, + "particle": "proton", + "RDD_model": "Gauss", # Options: Gauss, Geiss + "grid_size_um": 5.0, + "a0_nm": 8.0, + "use_beta": False, + "seed": None, # None for random seed + "SHOW_PLOT": False, + "PRINT_parameters": False, + "debug": False, + }, + "compare_initial_recombination": True, + "compare_continuous_beam": True, + "dose_rate_column": None, # None = auto-detect (will be null in YAML) + "energy_column": None, # None = auto-detect (will be null in YAML) + "ks_column": None, # None = auto-detect (will be null in YAML) + "plot_format": "png", + "plot_dpi": 300, + "plot_style": "whitegrid", + } + + with open(template_path, "w") as f: + yaml.dump(template, f, default_flow_style=False, sort_keys=False) diff --git a/validation/config.yaml b/validation/config.yaml new file mode 100644 index 0000000..22351cb --- /dev/null +++ b/validation/config.yaml @@ -0,0 +1,23 @@ +experimental_data_path: experimental-data/2019_DCPT/recombination_200V_data.csv +output_dir: results +simulation: + backend: cython + voltage_V: 200.0 + electrode_gap_cm: 0.2 + particle: proton + RDD_model: Gauss + grid_size_um: 5.0 + a0_nm: 8.0 + use_beta: false + seed: null + SHOW_PLOT: false + PRINT_parameters: false + debug: false +compare_initial_recombination: true +compare_continuous_beam: true +dose_rate_column: null +energy_column: null +ks_column: null +plot_format: png +plot_dpi: 300 +plot_style: whitegrid diff --git a/validation/config_example.yaml b/validation/config_example.yaml new file mode 100644 index 0000000..53c1d7a --- /dev/null +++ b/validation/config_example.yaml @@ -0,0 +1,57 @@ +# Example configuration file for IonTracks experimental data comparison +# +# Usage: +# 1. Copy this file to your desired location (e.g., validation/config.yaml or experimental-data/2019_DCPT/config.yaml) +# 2. Edit the paths and parameters as needed +# 3. Run: python -m validation.run_comparison validation/config.yaml +# +# Recommended locations: +# - validation/config.yaml - for project-wide configurations +# - experimental-data/DATASET_NAME/config.yaml - for experiment-specific configurations + +# Path to experimental data CSV file +# The CSV should contain columns for energy, dose rate, and k_s values +# See README.md for details on required column names +experimental_data_path: "experimental-data/2019_DCPT/recombination_200V_data.csv" + +# Output directory for results (CSV files and plots) +output_dir: "results" + +# Simulation configuration +simulation: + # Backend selection: cython (fastest), python, numba, cupy, parallel + backend: "cython" + + # Experimental setup parameters + voltage_V: 200.0 # Applied voltage in volts + electrode_gap_cm: 0.2 # Electrode gap in cm (typically 0.1-0.2 cm) + particle: "proton" # Particle type (proton, carbon, helium, etc.) + + # Physics model parameters + RDD_model: "Gauss" # Radial dose distribution model: "Gauss" or "Geiss" + grid_size_um: 5.0 # Spatial grid resolution in micrometers + # Smaller = more accurate but slower (typical: 3.0-5.0) + a0_nm: 8.0 # Track core radius parameter in nanometers + use_beta: false # Scale track core by β = v/c (usually false) + + # Computational options + seed: null # Random seed for reproducibility (null = random) + SHOW_PLOT: false # Show simulation plots during calculation + PRINT_parameters: false # Print detailed parameters + debug: false # Enable debug mode + +# Comparison options +compare_initial_recombination: true # Compare initial recombination (fast) +compare_continuous_beam: true # Compare general recombination (slow, may take hours) + +# Data column names (auto-detected if not specified) +# Leave as null to auto-detect, or specify exact column names +dose_rate_column: null # e.g., "dose_rate_air_Gy_s" or "dose_rate_water_Gy_s" +energy_column: null # e.g., "Energy_MeV" +ks_column: null # e.g., "k_s" + +# Plotting options +plot_format: "png" # Output format: "png", "pdf", "svg" +plot_dpi: 300 # Resolution for raster formats +plot_style: "whitegrid" # Matplotlib style + diff --git a/validation/data_utils.py b/validation/data_utils.py new file mode 100644 index 0000000..283a8e4 --- /dev/null +++ b/validation/data_utils.py @@ -0,0 +1,209 @@ +""" +Utilities for loading and validating experimental data. +""" + +from pathlib import Path +from typing import List, Optional, Tuple + +import pandas as pd +import numpy as np + + +REQUIRED_COLUMNS = { + "energy": ["Energy_MeV", "energy_MeV", "E_MeV", "energy"], + "ks": ["k_s", "ks", "recombination_factor", "collection_efficiency"], + "dose_rate_air": ["dose_rate_air_Gy_s", "dose_rate_air", "doserate_air_Gy_s"], + "dose_rate_water": ["dose_rate_water_Gy_s", "dose_rate_water", "doserate_water_Gy_s"], +} + + +def find_column(df: pd.DataFrame, possible_names: List[str]) -> Optional[str]: + """ + Find a column in DataFrame by trying multiple possible names. + + Parameters + ---------- + df : pd.DataFrame + DataFrame to search + possible_names : List[str] + List of possible column names + + Returns + ------- + Optional[str] + Found column name or None + """ + for name in possible_names: + if name in df.columns: + return name + return None + + +def validate_experimental_data( + df: pd.DataFrame, + energy_column: str, + ks_column: str, + dose_rate_column: str, +) -> tuple[bool, List[str]]: + """ + Validate experimental data DataFrame. + + Parameters + ---------- + df : pd.DataFrame + Experimental data DataFrame + energy_column : str + Name of energy column + ks_column : str + Name of k_s column + dose_rate_column : str + Name of dose rate column + + Returns + ------- + Tuple[bool, List[str]] + (is_valid, list_of_errors) + """ + errors = [] + + # Check required columns exist + required = [energy_column, ks_column, dose_rate_column] + missing = [col for col in required if col not in df.columns] + if missing: + errors.append(f"Missing required columns: {missing}") + return False, errors + + # Check for empty DataFrame + if len(df) == 0: + errors.append("DataFrame is empty") + return False, errors + + # Check for NaN values in required columns + for col in required: + nan_count = df[col].isna().sum() + if nan_count > 0: + errors.append(f"Column '{col}' contains {nan_count} NaN values") + + # Check data types and ranges + if not pd.api.types.is_numeric_dtype(df[energy_column]): + errors.append(f"Column '{energy_column}' must be numeric") + else: + if (df[energy_column] <= 0).any(): + errors.append(f"Column '{energy_column}' contains non-positive values") + + if not pd.api.types.is_numeric_dtype(df[ks_column]): + errors.append(f"Column '{ks_column}' must be numeric") + else: + if (df[ks_column] <= 0).any(): + errors.append(f"Column '{ks_column}' contains non-positive values") + + if not pd.api.types.is_numeric_dtype(df[dose_rate_column]): + errors.append(f"Column '{dose_rate_column}' must be numeric") + else: + if (df[dose_rate_column] < 0).any(): + errors.append(f"Column '{dose_rate_column}' contains negative values") + + return len(errors) == 0, errors + + +def load_experimental_data( + data_path: Path, + energy_column: Optional[str] = None, + ks_column: Optional[str] = None, + dose_rate_column: Optional[str] = None, +) -> pd.DataFrame: + """ + Load and validate experimental data from CSV file. + + Parameters + ---------- + data_path : Path + Path to CSV file with experimental data + energy_column : Optional[str] + Name of energy column (auto-detected if None) + ks_column : Optional[str] + Name of k_s column (auto-detected if None) + dose_rate_column : Optional[str] + Name of dose rate column (auto-detected if None) + + Returns + ------- + pd.DataFrame + Loaded and validated experimental data + + Raises + ------ + FileNotFoundError + If data file doesn't exist + ValueError + If data validation fails + """ + data_path = Path(data_path) + + if not data_path.exists(): + raise FileNotFoundError(f"Experimental data file not found: {data_path}") + + # Load CSV + try: + df = pd.read_csv(data_path) + except Exception as e: + raise ValueError(f"Failed to load CSV file: {e}") + + # Auto-detect columns if not specified + if energy_column is None: + energy_column = find_column(df, REQUIRED_COLUMNS["energy"]) + if energy_column is None: + raise ValueError( + f"Could not find energy column. Tried: {REQUIRED_COLUMNS['energy']}" + ) + + if ks_column is None: + ks_column = find_column(df, REQUIRED_COLUMNS["ks"]) + if ks_column is None: + raise ValueError( + f"Could not find k_s column. Tried: {REQUIRED_COLUMNS['ks']}" + ) + + if dose_rate_column is None: + # Try air first, then water + dose_rate_column = find_column(df, REQUIRED_COLUMNS["dose_rate_air"]) + if dose_rate_column is None: + dose_rate_column = find_column(df, REQUIRED_COLUMNS["dose_rate_water"]) + if dose_rate_column is None: + raise ValueError( + f"Could not find dose rate column. Tried: " + f"{REQUIRED_COLUMNS['dose_rate_air'] + REQUIRED_COLUMNS['dose_rate_water']}" + ) + + # Validate data + is_valid, errors = validate_experimental_data( + df, energy_column, ks_column, dose_rate_column + ) + + if not is_valid: + error_msg = "Data validation failed:\n" + "\n".join(f" - {e}" for e in errors) + raise ValueError(error_msg) + + # Standardize column names for internal use + df = df.copy() + df.rename( + columns={ + energy_column: "Energy_MeV", + ks_column: "k_s", + dose_rate_column: "dose_rate_Gy_s", + }, + inplace=True, + ) + + # Convert dose rate to Gy/min if needed (for continuous beam calculations) + if "dose_rate_Gy_min" not in df.columns: + df["dose_rate_Gy_min"] = df["dose_rate_Gy_s"] * 60 + + print(f"✓ Loaded {len(df)} experimental measurements") + print(f" Energy column: {energy_column} → Energy_MeV") + print(f" k_s column: {ks_column} → k_s") + print(f" Dose rate column: {dose_rate_column} → dose_rate_Gy_s") + print(f" Energies: {sorted(df['Energy_MeV'].unique())} MeV") + + return df + diff --git a/validation/plots.py b/validation/plots.py new file mode 100644 index 0000000..f128dd1 --- /dev/null +++ b/validation/plots.py @@ -0,0 +1,236 @@ +""" +Plotting functions for comparison results. +""" + +from pathlib import Path +from typing import Optional + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import seaborn as sns + +from validation.config import ComparisonConfig + + +def plot_initial_recombination( + comparison_df: pd.DataFrame, + output_path: Path, + config: ComparisonConfig, +): + """ + Plot initial recombination comparison. + + Parameters + ---------- + comparison_df : pd.DataFrame + Comparison DataFrame with initial recombination results + output_path : Path + Path to save the plot + config : ComparisonConfig + Configuration object + """ + if len(comparison_df) == 0: + print("Warning: No data to plot for initial recombination") + return + + sns.set_style(config.plot_style) + plt.rcParams["figure.dpi"] = config.plot_dpi + + fig, ax = plt.subplots(figsize=(10, 6)) + + ax.plot( + comparison_df["Energy_MeV"], + comparison_df["k_s_experimental"], + "o-", + label="Experimental", + markersize=10, + linewidth=2, + color="blue", + ) + ax.plot( + comparison_df["Energy_MeV"], + comparison_df["k_s_IonTracks"], + "s--", + label="IonTracks", + markersize=10, + linewidth=2, + color="red", + ) + + ax.set_xlabel("Proton Energy [MeV]", fontsize=12) + ax.set_ylabel("$k_s$", fontsize=12) + ax.set_title("Initial Recombination: Experiment vs. IonTracks", fontsize=14) + ax.legend(fontsize=11) + ax.grid(True, alpha=0.3) + + plt.tight_layout() + plt.savefig(output_path, format=config.plot_format, dpi=config.plot_dpi) + plt.close() + print(f"✓ Saved: {output_path.name}") + + +def plot_continuous_beam( + comparison_df: pd.DataFrame, + output_path: Path, + config: ComparisonConfig, +): + """ + Plot continuous beam comparison. + + Parameters + ---------- + comparison_df : pd.DataFrame + Comparison DataFrame with continuous beam results + output_path : Path + Path to save the plot + config : ComparisonConfig + Configuration object + """ + if len(comparison_df) == 0: + print("Warning: No data to plot for continuous beam") + return + + sns.set_style(config.plot_style) + plt.rcParams["figure.dpi"] = config.plot_dpi + + fig, ax = plt.subplots(figsize=(12, 7)) + + energies_MeV = sorted(comparison_df["Energy_MeV"].unique()) + colors = plt.cm.tab10(np.linspace(0, 1, len(energies_MeV))) + + for i, energy in enumerate(energies_MeV): + subset = comparison_df[comparison_df["Energy_MeV"] == energy] + subset = subset.sort_values("dose_rate_Gy_s") + + if len(subset) > 0: + ax.plot( + subset["dose_rate_Gy_s"], + subset["k_s"], + "o-", + label=f"Exp. {energy} MeV", + markersize=6, + linewidth=1.5, + color=colors[i], + alpha=0.7, + ) + ax.plot( + subset["dose_rate_Gy_s"], + subset["k_s_IonTracks"], + "s--", + label=f"IonTracks {energy} MeV", + markersize=6, + linewidth=1.5, + color=colors[i], + ) + + ax.set_xlabel("Dose Rate in Air [Gy/s]", fontsize=12) + ax.set_ylabel("$k_s$", fontsize=12) + ax.set_title("Continuous Beam: Experiment vs. IonTracks", fontsize=14) + ax.legend(ncol=2, fontsize=9) + ax.grid(True, alpha=0.3) + + plt.tight_layout() + plt.savefig(output_path, format=config.plot_format, dpi=config.plot_dpi) + plt.close() + print(f"✓ Saved: {output_path.name}") + + +def plot_relative_error( + comparison_df: pd.DataFrame, + output_path: Path, + config: ComparisonConfig, +): + """ + Plot relative error as a function of dose rate. + + Parameters + ---------- + comparison_df : pd.DataFrame + Comparison DataFrame with continuous beam results + output_path : Path + Path to save the plot + config : ComparisonConfig + Configuration object + """ + if len(comparison_df) == 0: + print("Warning: No data to plot for relative error") + return + + sns.set_style(config.plot_style) + plt.rcParams["figure.dpi"] = config.plot_dpi + + fig, ax = plt.subplots(figsize=(10, 6)) + + energies_MeV = sorted(comparison_df["Energy_MeV"].unique()) + colors = plt.cm.tab10(np.linspace(0, 1, len(energies_MeV))) + + for i, energy in enumerate(energies_MeV): + subset = comparison_df[comparison_df["Energy_MeV"] == energy] + subset = subset.sort_values("dose_rate_Gy_s") + + if len(subset) > 0: + ax.plot( + subset["dose_rate_Gy_s"], + subset["relative_error_%"], + "o-", + label=f"{energy} MeV", + markersize=6, + linewidth=1.5, + color=colors[i], + ) + + ax.axhline(y=0, color="r", linestyle=":", linewidth=2) + ax.set_xlabel("Dose Rate in Air [Gy/s]", fontsize=12) + ax.set_ylabel("Relative Error [%]", fontsize=12) + ax.set_title("IonTracks Relative Error", fontsize=14) + ax.legend(fontsize=11) + ax.grid(True, alpha=0.3) + + plt.tight_layout() + plt.savefig(output_path, format=config.plot_format, dpi=config.plot_dpi) + plt.close() + print(f"✓ Saved: {output_path.name}") + + +def generate_all_plots( + comparison_initial: Optional[pd.DataFrame], + comparison_continuous: Optional[pd.DataFrame], + config: ComparisonConfig, +): + """ + Generate all comparison plots. + + Parameters + ---------- + comparison_initial : Optional[pd.DataFrame] + Initial recombination comparison DataFrame + comparison_continuous : Optional[pd.DataFrame] + Continuous beam comparison DataFrame + config : ComparisonConfig + Configuration object + """ + print("\n" + "=" * 70) + print("GENERATING PLOTS") + print("=" * 70) + + if comparison_initial is not None and len(comparison_initial) > 0: + plot_initial_recombination( + comparison_initial, + config.output_dir / "initial_recombination_comparison.png", + config, + ) + + if comparison_continuous is not None and len(comparison_continuous) > 0: + plot_continuous_beam( + comparison_continuous, + config.output_dir / "continuous_beam_comparison.png", + config, + ) + + plot_relative_error( + comparison_continuous, + config.output_dir / "relative_error.png", + config, + ) + diff --git a/validation/run_comparison.py b/validation/run_comparison.py new file mode 100644 index 0000000..0dbf269 --- /dev/null +++ b/validation/run_comparison.py @@ -0,0 +1,90 @@ +""" +Main entry point for running experimental data comparisons. + +Usage: + python -m validation.run_comparison config.yaml +""" + +import argparse +import sys +from pathlib import Path + +from validation.config import ComparisonConfig, load_config, save_config_template +from validation.comparison import ComparisonRunner +from validation.plots import generate_all_plots + + +def main(): + """Main function.""" + parser = argparse.ArgumentParser( + description="Compare IonTracks simulations with experimental data" + ) + parser.add_argument( + "config", + type=str, + nargs="?", + help="Path to configuration YAML file (e.g., validation/config.yaml or experimental-data/2019_DCPT/config.yaml)", + ) + parser.add_argument( + "--create-template", + type=str, + metavar="PATH", + help="Create a template configuration file at the specified path", + ) + + args = parser.parse_args() + + # Create template if requested + if args.create_template: + template_path = Path(args.create_template) + save_config_template(template_path) + print(f"✓ Created template configuration at: {template_path}") + return + + # Load configuration + if not args.config: + parser.print_help() + print("\nError: Configuration file is required (or use --create-template)") + sys.exit(1) + + try: + config = load_config(args.config) + except Exception as e: + print(f"Error loading configuration: {e}") + sys.exit(1) + + # Run comparison + print("=" * 70) + print("IONTRACKS EXPERIMENTAL DATA COMPARISON") + print("=" * 70) + print(f"Experimental data: {config.experimental_data_path}") + print(f"Output directory: {config.output_dir}") + print(f"Backend: {config.simulation.backend}") + print(f"Voltage: {config.simulation.voltage_V} V") + print(f"Electrode gap: {config.simulation.electrode_gap_cm} cm") + print(f"Particle: {config.simulation.particle}") + + try: + runner = ComparisonRunner(config) + runner.run_all() + runner.save_results() + generate_all_plots( + runner.comparison_initial, + runner.comparison_continuous, + config, + ) + + print("\n" + "=" * 70) + print("COMPARISON COMPLETED SUCCESSFULLY") + print("=" * 70) + + except Exception as e: + print(f"\nError during comparison: {e}") + import traceback + traceback.print_exc() + sys.exit(1) + + +if __name__ == "__main__": + main() +