gammaloop_api/commands/

approach.rs

use gammalooprs::integrands::process::EvaluationTarget;
use std::{
    collections::{BTreeMap, BTreeSet},
    env, fs,
    io::{self, IsTerminal},
    path::{Path, PathBuf},
    sync::{
        atomic::{AtomicUsize, Ordering},
        Arc,
    },
};

use clap::Args;
use color_eyre::Result;
use colored::Colorize;
use eyre::{eyre, Context};
use gammalooprs::{
    integrands::{evaluation::EvaluationResult, process::ProcessIntegrand},
    model::Model,
    observables::events::{AdditionalWeightKey, GenericThresholdCountertermEventInfo},
    processes::ThresholdCountertermMetadataRegistry,
    utils::F,
};
use indicatif::ProgressBar;
use rayon::{prelude::*, ThreadPoolBuilder};
use schemars::JsonSchema;
use serde::{Deserialize, Serialize};
use spenso::algebra::complex::Complex;
use tracing::info;

use crate::{
    commands::evaluate_samples::{build_havana_sample, build_momentum_input},
    commands::CliArgumentMetadataExt,
    completion::CompletionArgExt,
    state::{ProcessRef, State},
    CLISettings,
};

const DEFAULT_MIN_ABS_T: f64 = 1.0e-6;
const DEFAULT_MAX_ABS_T: f64 = 1.0;

fn default_min_abs_t() -> f64 {
    DEFAULT_MIN_ABS_T
}

fn default_max_abs_t() -> f64 {
    DEFAULT_MAX_ABS_T
}

fn default_parameter_name() -> String {
    "t".to_string()
}

#[derive(Debug, Args, Serialize, Deserialize, Clone, JsonSchema, PartialEq)]
pub struct Approach {
    /// Process reference: `#<id>`, `name:<name>`, or `<id>/<name>`
    #[arg(
        short = 'p',
        long = "process",
        value_name = "PROCESS",
        completion_process_selector(crate::completion::SelectorKind::Any)
    )]
    pub process: Option<ProcessRef>,

    /// The integrand name to approach
    #[arg(
        short = 'i',
        long = "integrand-name",
        value_name = "NAME",
        completion_integrand_selector(crate::completion::SelectorKind::Any)
    )]
    pub integrand_name: Option<String>,

    /// Integration coordinates, or flattened loop momenta `(px, py, pz) ...` with `--momentum-space`
    #[arg(
        short = 'x',
        long = "point",
        num_args = 1..,
        value_name = "POINT",
        value_delimiter = ',',
        allow_negative_numbers = true,
    )]
    pub point: Vec<f64>,

    /// Direction vector in the same coordinate layout and dimension as `--point`
    #[arg(
        long = "approach-axis",
        value_name = "AXIS",
        num_args = 1..,
        allow_negative_numbers = true,
    )]
    pub approach_axes: Vec<String>,

    /// Human-readable label for each approach axis, in axis order
    #[arg(long = "axis-label", value_name = "LABEL")]
    #[serde(default)]
    pub axis_labels: Vec<String>,

    /// Number of points on each side of the midpoint
    #[arg(long = "n-points", value_name = "N")]
    pub n_points: usize,

    /// Do not evaluate the midpoint at t = 0
    #[arg(long)]
    #[serde(default)]
    pub skip_midpoint: bool,

    /// Use linear spacing in t
    #[arg(long, conflicts_with = "logarithmic")]
    pub linear: bool,

    /// Use logarithmic spacing in |t|
    #[arg(long, conflicts_with = "linear")]
    pub logarithmic: bool,

    /// Smallest non-zero |t| for logarithmic spacing
    #[arg(
        long = "min-abs-t",
        default_value_t = DEFAULT_MIN_ABS_T,
        allow_negative_numbers = true
    )]
    #[serde(default = "default_min_abs_t")]
    pub min_abs_t: f64,

    /// Largest |t| evaluated on either side of the midpoint
    #[arg(
        long = "max-abs-t",
        default_value_t = DEFAULT_MAX_ABS_T,
        allow_negative_numbers = true
    )]
    #[serde(default = "default_max_abs_t")]
    pub max_abs_t: f64,

    /// Name of the physical parameter t in point(t) = midpoint + t * axis
    #[arg(long = "parameter-name", default_value = "t", value_name = "NAME")]
    #[serde(default = "default_parameter_name")]
    pub parameter_name: String,

    /// Number of worker threads for approach evaluations
    #[arg(long = "n-cores", value_name = "N")]
    pub n_cores: Option<usize>,

    /// Force arbitrary-precision (Arb) internal evaluation; results remain f64
    #[arg(short = 'f', long = "use_arb_prec")]
    pub use_arb_prec: bool,

    /// Interpret `--point` and every axis as spatial loop-momentum triplets
    #[arg(short = 'm', long)]
    pub momentum_space: bool,

    /// The discrete dimensions of the sample
    #[arg(
        short = 'd',
        long = "discrete-dim",
        value_name = "DIMS",
        num_args = 1..,
        value_delimiter = ',',
        conflicts_with_all = ["graph_id", "orientation_id"],
    )]
    pub discrete_dim: Vec<usize>,

    /// Select a specific graph in momentum-space approach
    #[arg(long = "graph-id", value_name = "GRAPH_ID")]
    pub graph_id: Option<usize>,

    /// Select a specific orientation of the selected graph in momentum-space approach
    #[arg(
        long = "orientation-id",
        value_name = "ORIENTATION_ID",
        cli_requires("graph_id"),
        conflicts_with = "discrete_dim"
    )]
    pub orientation_id: Option<usize>,

    /// Write sampled approach points and run metadata to this JSON file
    #[arg(long = "output-results", value_hint = clap::ValueHint::FilePath)]
    pub output_results: Option<PathBuf>,
}

impl Default for Approach {
    fn default() -> Self {
        Self {
            process: None,
            integrand_name: None,
            point: Vec::new(),
            approach_axes: Vec::new(),
            axis_labels: Vec::new(),
            n_points: 0,
            skip_midpoint: false,
            linear: false,
            logarithmic: false,
            min_abs_t: DEFAULT_MIN_ABS_T,
            max_abs_t: DEFAULT_MAX_ABS_T,
            parameter_name: default_parameter_name(),
            n_cores: None,
            use_arb_prec: false,
            momentum_space: false,
            discrete_dim: Vec::new(),
            graph_id: None,
            orientation_id: None,
            output_results: None,
        }
    }
}

#[derive(Debug, Clone, Copy, Serialize, Deserialize, PartialEq, Eq)]
#[serde(rename_all = "snake_case")]
enum ApproachSpacing {
    Linear,
    Logarithmic,
}

#[derive(Debug, Clone)]
struct ApproachJob {
    index: usize,
    axis_index: usize,
    axis_point_index: usize,
    t: f64,
    point: Vec<f64>,
    skip_reason: Option<String>,
}

#[derive(Debug, Serialize)]
struct ApproachOutput {
    schema_version: u32,
    command: ApproachCommandJson,
    process: ApproachProcessJson,
    integrand: ApproachIntegrandJson,
    space: String,
    base_point: Vec<f64>,
    parameter_name: String,
    axes: Vec<ApproachAxisJson>,
    spacing: ApproachSpacingJson,
    n_cores: usize,
    points_per_axis: usize,
    evaluated_points: usize,
    skipped_points: usize,
    points: Vec<ApproachPointRecord>,
}

#[derive(Debug, Serialize)]
struct ApproachCommandJson {
    name: &'static str,
    use_arb_prec: bool,
    skip_midpoint: bool,
    graph_id: Option<usize>,
    graph_name: Option<String>,
    orientation_id: Option<usize>,
    discrete_dim: Vec<usize>,
}

#[derive(Debug, Serialize)]
struct ApproachProcessJson {
    id: usize,
    name: String,
}

#[derive(Debug, Serialize)]
struct ApproachIntegrandJson {
    name: String,
    kind: String,
    #[serde(skip_serializing_if = "Vec::is_empty")]
    threshold_counterterms: Vec<ApproachThresholdCountertermRegistryJson>,
}

#[derive(Debug, Serialize)]
struct ApproachThresholdCountertermRegistryJson {
    graph_id: usize,
    registry: ThresholdCountertermMetadataRegistry,
}

#[derive(Debug, Serialize)]
struct ApproachSpacingJson {
    kind: ApproachSpacing,
    n_points: usize,
    min_abs_t: Option<f64>,
    max_abs_t: f64,
    t_values: Vec<f64>,
}

#[derive(Debug, Clone, Serialize)]
struct ApproachAxisJson {
    label: String,
    vector: Vec<f64>,
}

#[derive(Debug, Serialize)]
struct ApproachPointRecord {
    index: usize,
    axis_index: usize,
    axis_point_index: usize,
    t: f64,
    point: Vec<f64>,
    status: &'static str,
    #[serde(skip_serializing_if = "Option::is_none")]
    skip_reason: Option<String>,
    #[serde(skip_serializing_if = "Option::is_none")]
    evaluation: Option<ApproachEvaluationRecord>,
}

#[derive(Debug, Serialize)]
struct ApproachEvaluationRecord {
    integrand_result: ComplexJson,
    parameterization_jacobian: Option<f64>,
    integrator_weight: f64,
    total_weight: ComplexJson,
    event_weight_sum: ComplexJson,
    additional_contribution_sums: BTreeMap<String, ComplexJson>,
    #[serde(skip_serializing_if = "Option::is_none")]
    threshold_counterterm_summary: Option<ApproachThresholdCountertermSummaryJson>,
    contributions: Vec<ContributionRecord>,
    events: Vec<EventRecord>,
    metadata: EvaluationMetadataJson,
}

#[derive(Debug, Serialize)]
struct EvaluationMetadataJson {
    generated_event_count: usize,
    accepted_event_count: usize,
    is_nan: bool,
    total_time_seconds: f64,
    parameterization_time_seconds: f64,
    integrand_evaluation_time_seconds: f64,
    evaluator_evaluation_time_seconds: f64,
    event_processing_time_seconds: f64,
}

#[derive(Debug, Serialize)]
struct EventRecord {
    event_group_index: usize,
    event_index: usize,
    graph_id: usize,
    graph_name: Option<String>,
    graph_group_id: Option<usize>,
    orientation_id: Option<usize>,
    cut_id: usize,
    cut_edges: Vec<usize>,
    sampling_channel_id: Option<usize>,
    lmb_sample_id: Option<usize>,
    weight: ComplexJson,
    #[serde(skip_serializing_if = "Option::is_none")]
    full_multiplicative_factor: Option<ComplexJson>,
    display_only_weights: BTreeMap<String, ComplexJson>,
    #[serde(skip_serializing_if = "Option::is_none")]
    threshold_counterterms: Option<GenericThresholdCountertermEventInfo<f64>>,
}

#[derive(Debug, Serialize)]
struct ApproachThresholdCountertermSummaryJson {
    original: ComplexJson,
    counterterm_sum: ComplexJson,
    reconstructed_total: ComplexJson,
    decomposed_event_sum: ComplexJson,
    closure: ComplexJson,
    components: Vec<ApproachThresholdCountertermComponentSumJson>,
}

#[derive(Debug, Serialize)]
struct ApproachThresholdCountertermComponentSumJson {
    graph_id: usize,
    component_id: usize,
    weighted: ComplexJson,
    occurrence_count: usize,
    skipped_count: usize,
}

#[derive(Debug)]
struct ThresholdCountertermComponentAccumulator {
    weighted: Complex<F<f64>>,
    occurrence_count: usize,
    skipped_count: usize,
}

#[derive(Debug)]
struct ThresholdCountertermSummaryAccumulator {
    decomposed_event_count: usize,
    original: Complex<F<f64>>,
    decomposed_event_sum: Complex<F<f64>>,
    components: BTreeMap<(usize, usize), ThresholdCountertermComponentAccumulator>,
}

impl Default for ThresholdCountertermSummaryAccumulator {
    fn default() -> Self {
        Self {
            decomposed_event_count: 0,
            original: zero_complex(),
            decomposed_event_sum: zero_complex(),
            components: BTreeMap::new(),
        }
    }
}

impl ThresholdCountertermSummaryAccumulator {
    fn add(
        &mut self,
        graph_id: usize,
        event_weight: &Complex<F<f64>>,
        decomposition: &GenericThresholdCountertermEventInfo<f64>,
    ) {
        self.decomposed_event_count += 1;
        self.decomposed_event_sum += *event_weight;
        self.original += decomposition.original;
        for component in &decomposition.components {
            let accumulator = self
                .components
                .entry((graph_id, component.component_id))
                .or_insert_with(|| ThresholdCountertermComponentAccumulator {
                    weighted: zero_complex(),
                    occurrence_count: 0,
                    skipped_count: 0,
                });
            accumulator.weighted += component.weighted;
            accumulator.occurrence_count += 1;
            accumulator.skipped_count += usize::from(component.evaluation_skipped);
        }
    }

    fn into_json(self) -> Option<ApproachThresholdCountertermSummaryJson> {
        if self.decomposed_event_count == 0 {
            return None;
        }
        let counterterm_sum = self
            .components
            .values()
            .fold(zero_complex(), |total, component| {
                total + component.weighted
            });
        let reconstructed_total = self.original + counterterm_sum;
        let closure = self.decomposed_event_sum - reconstructed_total;
        Some(ApproachThresholdCountertermSummaryJson {
            original: complex_json(self.original),
            counterterm_sum: complex_json(counterterm_sum),
            reconstructed_total: complex_json(reconstructed_total),
            decomposed_event_sum: complex_json(self.decomposed_event_sum),
            closure: complex_json(closure),
            components: self
                .components
                .into_iter()
                .map(|((graph_id, component_id), component)| {
                    ApproachThresholdCountertermComponentSumJson {
                        graph_id,
                        component_id,
                        weighted: complex_json(component.weighted),
                        occurrence_count: component.occurrence_count,
                        skipped_count: component.skipped_count,
                    }
                })
                .collect(),
        })
    }
}

#[derive(Debug, Serialize)]
struct ContributionRecord {
    label: String,
    contribution: String,
    graph_id: usize,
    graph_name: Option<String>,
    graph_group_id: Option<usize>,
    orientation_id: Option<usize>,
    cut_id: usize,
    cut_edges: Vec<usize>,
    lmb_sample_id: Option<usize>,
    weight: ComplexJson,
}

#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord)]
struct ContributionKey {
    contribution: String,
    graph_id: usize,
    graph_name: Option<String>,
    graph_group_id: Option<usize>,
    orientation_id: Option<usize>,
    cut_id: usize,
    cut_edges: Vec<usize>,
    lmb_sample_id: Option<usize>,
}

#[derive(Debug, Clone, Copy, Serialize)]
struct ComplexJson {
    re: f64,
    im: f64,
}

impl Approach {
    pub fn run(&self, state: &mut State, cli_settings: &CLISettings) -> Result<PathBuf> {
        let output_path = self
            .output_results
            .clone()
            .unwrap_or_else(|| PathBuf::from("approach_result.json"));
        cli_settings
            .ensure_write_target_outside_active_state(&output_path, "write approach results")?;
        self.validate_cli()?;
        let (process_id, integrand_name) =
            state.find_integrand_ref(self.process.as_ref(), self.integrand_name.as_ref())?;
        let process_name = state.process_list.processes[process_id]
            .definition
            .folder_name
            .clone();
        let model = state.resolve_model_for_integrand(process_id, &integrand_name)?;
        let base_integrand = state
            .process_list
            .get_integrand(process_id, &integrand_name)?
            .require_generated()?
            .clone();

        if !self.momentum_space && (self.graph_id.is_some() || self.orientation_id.is_some()) {
            return Err(eyre!(
                "Graph and orientation selectors are only supported in momentum-space approach."
            ));
        }
        if !self.momentum_space {
            let expected_dimension =
                base_integrand.expected_x_space_dimension(&self.discrete_dim)?;
            if self.point.len() != expected_dimension {
                return Err(eyre!(
                    "Expected {} x-space coordinates for this integrand selection, got {}.",
                    expected_dimension,
                    self.point.len()
                ));
            }
        }

        let graph_name = self.resolve_graph_name(&base_integrand, &integrand_name)?;
        let axes = self.parse_axes()?;
        let spacing = self.spacing()?;
        let t_values = self.t_values(spacing)?;
        let points_per_axis = t_values.len();
        let jobs = self.build_jobs(&axes, &t_values);
        let n_cores = self
            .n_cores
            .unwrap_or(cli_settings.global.n_cores.integrate);
        if n_cores == 0 {
            return Err(eyre!("--n-cores must be at least 1."));
        }

        let total_jobs = jobs.len();
        let progress = approach_progress_bar(total_jobs as u64);
        let pool = ThreadPoolBuilder::new()
            .num_threads(n_cores)
            .build()
            .context("failed to build approach evaluation thread pool")?;
        let worker_count = n_cores.min(total_jobs.max(1));
        let worker_inputs = (0..worker_count)
            .map(|worker_index| {
                let worker_jobs = jobs
                    .iter()
                    .skip(worker_index)
                    .step_by(worker_count)
                    .cloned()
                    .collect::<Vec<_>>();
                let mut integrand = base_integrand.clone();
                force_event_output(&mut integrand);
                (integrand, model.clone(), worker_jobs)
            })
            .collect::<Vec<_>>();
        let completed_jobs = Arc::new(AtomicUsize::new(0));
        let completed_workers = Arc::new(AtomicUsize::new(0));
        let results = pool.install(|| {
            worker_inputs
                .into_par_iter()
                .map(
                    |(mut integrand, model, worker_jobs)| -> Result<Vec<ApproachPointRecord>> {
                        integrand
                            .warm_up(&model)
                            .context("failed to warm up approach worker integrand")?;
                        let mut records = Vec::with_capacity(worker_jobs.len());
                        for job in worker_jobs {
                            let record = self.evaluate_job(
                                &job,
                                &mut integrand,
                                &model,
                                graph_name.as_deref(),
                                axes.len(),
                                t_values.len(),
                            );
                            let completed_total =
                                completed_jobs.fetch_add(1, Ordering::Relaxed) + 1;
                            progress.inc(1);
                            progress.set_message(progress_message(
                                job.axis_index,
                                axes.len(),
                                job.axis_point_index,
                                t_values.len(),
                                completed_total,
                                total_jobs,
                            ));
                            records.push(record?);
                        }
                        let completed_worker_count =
                            completed_workers.fetch_add(1, Ordering::Relaxed) + 1;
                        progress.set_message(format!(
                            "{} {}",
                            "collecting worker results".bright_magenta().bold(),
                            format!("{completed_worker_count:>2}/{worker_count:<2}")
                                .bright_white()
                                .bold()
                        ));
                        if completed_worker_count == worker_count {
                            progress.println(
                                "Approach evaluations complete; collecting worker results"
                                    .bright_blue()
                                    .bold()
                                    .to_string(),
                            );
                        }
                        Ok(records)
                    },
                )
                .collect::<Result<Vec<_>>>()
        });
        progress.finish_with_message(
            "approach evaluations complete; finalizing output"
                .bright_green()
                .bold()
                .to_string(),
        );

        info!(
            "{} {}",
            "Finalizing approach results".bright_blue().bold(),
            "(sorting records and preparing JSON)".bright_black()
        );
        let mut point_records = results?
            .into_iter()
            .flatten()
            .collect::<Vec<ApproachPointRecord>>();
        point_records.sort_by_key(|record| record.index);

        let evaluated_points = point_records
            .iter()
            .filter(|record| record.status == "evaluated")
            .count();
        let skipped_points = point_records.len() - evaluated_points;
        let output = ApproachOutput {
            schema_version: 3,
            command: ApproachCommandJson {
                name: "approach",
                use_arb_prec: self.use_arb_prec,
                skip_midpoint: self.skip_midpoint,
                graph_id: self.graph_id,
                graph_name: graph_name.clone(),
                orientation_id: self.orientation_id,
                discrete_dim: self.discrete_dim.clone(),
            },
            process: ApproachProcessJson {
                id: process_id,
                name: process_name,
            },
            integrand: ApproachIntegrandJson {
                name: integrand_name,
                kind: base_integrand.kind_name().to_string(),
                threshold_counterterms: threshold_counterterm_registries(
                    &base_integrand,
                    self.graph_id,
                ),
            },
            space: if self.momentum_space {
                "momentum".to_string()
            } else {
                "coordinate".to_string()
            },
            base_point: self.point.clone(),
            parameter_name: self.parameter_name.trim().to_string(),
            axes,
            spacing: ApproachSpacingJson {
                kind: spacing,
                n_points: self.n_points,
                min_abs_t: matches!(spacing, ApproachSpacing::Logarithmic)
                    .then_some(self.min_abs_t),
                max_abs_t: self.max_abs_t,
                t_values,
            },
            n_cores,
            points_per_axis,
            evaluated_points,
            skipped_points,
            points: point_records,
        };

        let relative_output = relative_path_display(&output_path);
        info!(
            "{} {}",
            "Writing approach JSON to".bright_blue().bold(),
            relative_output.bright_cyan()
        );
        if let Some(parent) = output_path.parent() {
            if !parent.as_os_str().is_empty() {
                fs::create_dir_all(parent).with_context(|| {
                    format!(
                        "failed to create approach output directory '{}'",
                        parent.display()
                    )
                })?;
            }
        }
        let file = fs::File::create(&output_path).with_context(|| {
            format!(
                "failed to create approach output file '{}'",
                output_path.display()
            )
        })?;
        serde_json::to_writer_pretty(file, &output).with_context(|| {
            format!(
                "failed to serialize approach output to '{}'",
                output_path.display()
            )
        })?;

        let pdf_path = output_path.with_extension("pdf");
        let relative_pdf = relative_path_display(&pdf_path);
        let plot_command = format!(
            "python3 assets/plot_approach_result.py {} --output {}",
            shell_quote(&relative_output),
            shell_quote(&relative_pdf)
        );
        info!(
            "{} {}",
            "Approach results written to".bright_green().bold(),
            relative_output.bright_cyan()
        );
        info!(
            "{} {}",
            "Plot with:".bright_green().bold(),
            plot_command.bright_cyan()
        );

        Ok(output_path)
    }

    fn validate_cli(&self) -> Result<()> {
        if self.point.is_empty() {
            return Err(eyre!("approach requires a midpoint supplied with --point."));
        }
        if let Some((index, value)) = self
            .point
            .iter()
            .enumerate()
            .find(|(_, value)| !value.is_finite())
        {
            return Err(eyre!(
                "Approach midpoint component {index} must be finite, got {value}."
            ));
        }
        if self.approach_axes.is_empty() {
            return Err(eyre!(
                "approach requires at least one --approach-axis value."
            ));
        }
        if !self.axis_labels.is_empty() && self.axis_labels.len() != self.approach_axes.len() {
            return Err(eyre!(
                "Expected either no --axis-label values or one per approach axis ({}), got {}.",
                self.approach_axes.len(),
                self.axis_labels.len(),
            ));
        }
        let mut labels = BTreeSet::new();
        for label in &self.axis_labels {
            let label = label.trim();
            if label.is_empty() {
                return Err(eyre!("Approach axis labels must not be blank."));
            }
            if !labels.insert(label) {
                return Err(eyre!("Approach axis label '{label}' is duplicated."));
            }
        }
        if self.parameter_name.trim().is_empty() {
            return Err(eyre!("--parameter-name must not be blank."));
        }
        if self.n_points == 0 {
            return Err(eyre!("--n-points must be at least 1."));
        }
        if !self.max_abs_t.is_finite() || self.max_abs_t <= 0.0 {
            return Err(eyre!("--max-abs-t must be finite and greater than 0."));
        }
        if !self.min_abs_t.is_finite() || self.min_abs_t <= 0.0 {
            return Err(eyre!("--min-abs-t must be finite and greater than 0."));
        }
        if self.logarithmic && self.min_abs_t > self.max_abs_t {
            return Err(eyre!(
                "--min-abs-t must be at most --max-abs-t for logarithmic spacing."
            ));
        }
        if self.n_cores == Some(0) {
            return Err(eyre!("--n-cores must be at least 1."));
        }
        Ok(())
    }

    fn spacing(&self) -> Result<ApproachSpacing> {
        if self.linear && self.logarithmic {
            return Err(eyre!("--linear and --logarithmic are mutually exclusive."));
        }
        Ok(if self.logarithmic {
            ApproachSpacing::Logarithmic
        } else {
            ApproachSpacing::Linear
        })
    }

    fn parse_axes(&self) -> Result<Vec<ApproachAxisJson>> {
        self.approach_axes
            .iter()
            .enumerate()
            .map(|(axis_index, raw_axis)| {
                Ok(ApproachAxisJson {
                    label: self
                        .axis_labels
                        .get(axis_index)
                        .map(|label| label.trim().to_string())
                        .unwrap_or_else(|| format!("axis {}", axis_index + 1)),
                    vector: parse_axis(raw_axis, self.point.len())?,
                })
            })
            .collect()
    }

    fn resolve_graph_name(
        &self,
        integrand: &ProcessIntegrand,
        integrand_name: &str,
    ) -> Result<Option<String>> {
        let Some(graph_id) = self.graph_id else {
            return Ok(None);
        };
        let graph_name = integrand.graph_name_by_id(graph_id).ok_or_else(|| {
            eyre!(
                "Graph id {} is out of range for integrand '{}'; it has {} graphs.",
                graph_id,
                integrand_name,
                integrand.graph_count()
            )
        })?;
        Ok(Some(graph_name.to_string()))
    }

    fn t_values(&self, spacing: ApproachSpacing) -> Result<Vec<f64>> {
        let magnitudes = match spacing {
            ApproachSpacing::Linear => (1..=self.n_points)
                .map(|index| self.max_abs_t * index as f64 / self.n_points as f64)
                .collect::<Vec<_>>(),
            ApproachSpacing::Logarithmic => {
                if self.n_points == 1 {
                    vec![self.max_abs_t]
                } else {
                    let log_min = self.min_abs_t.ln();
                    let log_max = self.max_abs_t.ln();
                    let denom = (self.n_points - 1) as f64;
                    (0..self.n_points)
                        .map(|index| (log_min + (log_max - log_min) * index as f64 / denom).exp())
                        .collect()
                }
            }
        };
        let mut t_values = magnitudes
            .iter()
            .rev()
            .map(|value| -*value)
            .collect::<Vec<_>>();
        if !self.skip_midpoint {
            t_values.push(0.0);
        }
        t_values.extend(magnitudes);
        Ok(t_values)
    }

    fn build_jobs(&self, axes: &[ApproachAxisJson], t_values: &[f64]) -> Vec<ApproachJob> {
        let mut jobs = Vec::with_capacity(axes.len() * t_values.len());
        for (axis_index, axis) in axes.iter().enumerate() {
            for (axis_point_index, t) in t_values.iter().copied().enumerate() {
                let point = self
                    .point
                    .iter()
                    .zip(&axis.vector)
                    .map(|(center, direction)| center + t * direction)
                    .collect::<Vec<_>>();
                let skip_reason = (!self.momentum_space)
                    .then(|| coordinate_skip_reason(&point))
                    .flatten();
                jobs.push(ApproachJob {
                    index: jobs.len(),
                    axis_index,
                    axis_point_index,
                    t,
                    point,
                    skip_reason,
                });
            }
        }
        jobs
    }

    fn evaluate_job(
        &self,
        job: &ApproachJob,
        integrand: &mut ProcessIntegrand,
        model: &Model,
        graph_name: Option<&str>,
        _axis_count: usize,
        _axis_point_count: usize,
    ) -> Result<ApproachPointRecord> {
        if let Some(skip_reason) = &job.skip_reason {
            return Ok(ApproachPointRecord {
                index: job.index,
                axis_index: job.axis_index,
                axis_point_index: job.axis_point_index,
                t: job.t,
                point: job.point.clone(),
                status: "skipped",
                skip_reason: Some(skip_reason.clone()),
                evaluation: None,
            });
        }

        let mut samples = if self.momentum_space {
            let input = build_momentum_input(
                integrand,
                &job.point,
                1.0,
                &self.discrete_dim,
                graph_name,
                self.orientation_id,
            )?;
            integrand
                .evaluate_momentum_configurations_raw(model, &[input], self.use_arb_prec)?
                .samples
        } else {
            let sample = build_havana_sample(integrand, &job.point, &self.discrete_dim, 1.0)?;
            integrand
                .evaluate_samples_raw(
                    EvaluationTarget::Physical(model),
                    &[sample],
                    1,
                    self.use_arb_prec,
                    false,
                    Complex::new_zero(),
                )?
                .samples
        };
        let evaluation = samples
            .pop()
            .ok_or_else(|| eyre!("approach evaluation did not return a sample"))?;
        Ok(ApproachPointRecord {
            index: job.index,
            axis_index: job.axis_index,
            axis_point_index: job.axis_point_index,
            t: job.t,
            point: job.point.clone(),
            status: "evaluated",
            skip_reason: None,
            evaluation: Some(approach_evaluation_record(integrand, evaluation)?),
        })
    }
}

fn parse_axis(raw_axis: &str, expected_dimension: usize) -> Result<Vec<f64>> {
    let trimmed = raw_axis
        .trim()
        .trim_start_matches('[')
        .trim_end_matches(']');
    let axis = trimmed
        .split(',')
        .map(str::trim)
        .filter(|component| !component.is_empty())
        .map(|component| {
            component.parse::<f64>().map_err(|err| {
                eyre!(
                    "Could not parse approach-axis component '{}' as a float: {}",
                    component,
                    err
                )
            })
        })
        .collect::<Result<Vec<_>>>()?;
    if axis.len() != expected_dimension {
        return Err(eyre!(
            "Approach axis has {} components, but the point has {} components.",
            axis.len(),
            expected_dimension
        ));
    }
    if let Some((index, value)) = axis
        .iter()
        .enumerate()
        .find(|(_, value)| !value.is_finite())
    {
        return Err(eyre!(
            "Approach axis component {index} must be finite, got {value}."
        ));
    }
    if axis.iter().all(|component| *component == 0.0) {
        return Err(eyre!("Approach axis must not be the zero vector."));
    }
    Ok(axis)
}

fn coordinate_skip_reason(point: &[f64]) -> Option<String> {
    let invalid = point
        .iter()
        .enumerate()
        .filter(|(_, value)| !(0.0..=1.0).contains(*value))
        .map(|(index, value)| format!("{index}:{value:+.16e}"))
        .collect::<Vec<_>>();
    (!invalid.is_empty()).then(|| {
        format!(
            "coordinate-space point lies outside the unit hypercube at component(s) {}",
            invalid.join(", ")
        )
    })
}

fn force_event_output(integrand: &mut ProcessIntegrand) {
    let settings = integrand.get_mut_settings();
    settings.general.generate_events = true;
    settings.general.store_additional_weights_in_event = true;
}

fn threshold_counterterm_registries(
    integrand: &ProcessIntegrand,
    selected_graph_id: Option<usize>,
) -> Vec<ApproachThresholdCountertermRegistryJson> {
    let collect = |registries: Vec<Option<&ThresholdCountertermMetadataRegistry>>| {
        registries
            .into_iter()
            .enumerate()
            .filter(|(graph_id, _)| selected_graph_id.is_none_or(|selected| selected == *graph_id))
            .filter_map(|(graph_id, registry)| {
                registry
                    .cloned()
                    .map(|registry| ApproachThresholdCountertermRegistryJson { graph_id, registry })
            })
            .collect()
    };
    match integrand {
        ProcessIntegrand::Amplitude(integrand) => collect(
            integrand
                .data
                .graph_terms
                .iter()
                .map(|graph| graph.threshold_counterterm_metadata())
                .collect(),
        ),
        ProcessIntegrand::CrossSection(integrand) => collect(
            integrand
                .data
                .graph_terms
                .iter()
                .map(|graph| graph.threshold_counterterm_metadata())
                .collect(),
        ),
    }
}

fn approach_evaluation_record(
    integrand: &ProcessIntegrand,
    evaluation: EvaluationResult,
) -> Result<ApproachEvaluationRecord> {
    let parameterization_jacobian = evaluation.parameterization_jacobian.map(|value| value.0);
    let integrator_weight = evaluation.integrator_weight.0;
    let total_scale = parameterization_jacobian.unwrap_or(1.0) * integrator_weight;
    let total_weight = Complex::new(
        F(evaluation.integrand_result.re.0 * total_scale),
        F(evaluation.integrand_result.im.0 * total_scale),
    );

    let mut event_weight_sum = zero_complex();
    let mut additional_contribution_sums = BTreeMap::<String, Complex<F<f64>>>::new();
    let mut contribution_sums = BTreeMap::<ContributionKey, Complex<F<f64>>>::new();
    let mut threshold_counterterm_summary = ThresholdCountertermSummaryAccumulator::default();
    let mut events = Vec::new();
    let parameterization_settings = integrand
        .get_settings()
        .sampling
        .get_parameterization_settings()
        .unwrap_or_default();

    for (event_group_index, event_group) in evaluation.event_groups.iter().enumerate() {
        for (event_index, event) in event_group.iter().enumerate() {
            event_weight_sum += event.weight;
            let graph_id = event.cut_info.graph_id;
            let graph_name = integrand
                .graph_name_by_id(graph_id)
                .map(ToString::to_string);
            let graph_group_id = integrand.graph_group_id_by_graph_id(graph_id);
            let cut_id = event.cut_info.cut_id;
            let cut_edges = integrand.cut_edge_ids(graph_id, cut_id).unwrap_or_default();
            let sampling_channel_id = event.cut_info.sampling_channel_id;
            let lmb_sample_id = sampling_channel_id
                .map(|channel_id| {
                    integrand.lmb_sample_id_for_channel(
                        graph_id,
                        channel_id,
                        &parameterization_settings,
                    )
                })
                .transpose()?
                .flatten();

            let event_key = ContributionKey {
                contribution: "event_weight".to_string(),
                graph_id,
                graph_name: graph_name.clone(),
                graph_group_id,
                orientation_id: event.cut_info.orientation_id,
                cut_id,
                cut_edges: cut_edges.clone(),
                lmb_sample_id,
            };
            add_contribution(&mut contribution_sums, event_key, event.weight);

            let full_multiplicative_factor = event
                .additional_weights
                .weights
                .get(&AdditionalWeightKey::FullMultiplicativeFactor)
                .cloned();
            let has_display_only_weights = event
                .additional_weights
                .weights
                .keys()
                .any(|key| *key != AdditionalWeightKey::FullMultiplicativeFactor);
            if has_display_only_weights && full_multiplicative_factor.is_none() {
                return Err(eyre!(
                    "Approach event for graph {} cut {} has display-only weights but no full multiplicative factor.",
                    graph_id,
                    cut_id,
                ));
            }

            let mut display_only_weights = BTreeMap::new();
            for (key, value) in &event.additional_weights.weights {
                let Some(label) = additional_weight_key_label(*key) else {
                    continue;
                };
                let normalized = *value
                    * *full_multiplicative_factor
                        .as_ref()
                        .expect("display-only weights require a full multiplicative factor");
                *additional_contribution_sums
                    .entry(label.clone())
                    .or_insert_with(zero_complex) += normalized;
                display_only_weights.insert(label.clone(), complex_json(*value));
                let contribution_key = ContributionKey {
                    contribution: label,
                    graph_id,
                    graph_name: graph_name.clone(),
                    graph_group_id,
                    orientation_id: event.cut_info.orientation_id,
                    cut_id,
                    cut_edges: cut_edges.clone(),
                    lmb_sample_id,
                };
                add_contribution(&mut contribution_sums, contribution_key, normalized);
            }

            if let Some(decomposition) = &event.additional_weights.threshold_counterterms {
                threshold_counterterm_summary.add(graph_id, &event.weight, decomposition);
            }

            events.push(EventRecord {
                event_group_index,
                event_index,
                graph_id,
                graph_name,
                graph_group_id,
                orientation_id: event.cut_info.orientation_id,
                cut_id,
                cut_edges,
                sampling_channel_id,
                lmb_sample_id,
                weight: complex_json(event.weight),
                full_multiplicative_factor: full_multiplicative_factor.map(complex_json),
                display_only_weights,
                threshold_counterterms: event.additional_weights.threshold_counterterms.clone(),
            });
        }
    }
    let threshold_counterterm_summary = threshold_counterterm_summary.into_json();

    Ok(ApproachEvaluationRecord {
        integrand_result: complex_json(evaluation.integrand_result),
        parameterization_jacobian,
        integrator_weight,
        total_weight: complex_json(total_weight),
        event_weight_sum: complex_json(event_weight_sum),
        additional_contribution_sums: additional_contribution_sums
            .into_iter()
            .map(|(key, value)| (key, complex_json(value)))
            .collect(),
        threshold_counterterm_summary,
        contributions: contribution_sums
            .into_iter()
            .map(|(key, weight)| contribution_record(key, weight))
            .collect(),
        events,
        metadata: EvaluationMetadataJson {
            generated_event_count: evaluation.evaluation_metadata.generated_event_count,
            accepted_event_count: evaluation.evaluation_metadata.accepted_event_count,
            is_nan: evaluation.evaluation_metadata.is_nan,
            total_time_seconds: evaluation.evaluation_metadata.total_timing.as_secs_f64(),
            parameterization_time_seconds: evaluation
                .evaluation_metadata
                .parameterization_time
                .as_secs_f64(),
            integrand_evaluation_time_seconds: evaluation
                .evaluation_metadata
                .integrand_evaluation_time
                .as_secs_f64(),
            evaluator_evaluation_time_seconds: evaluation
                .evaluation_metadata
                .evaluator_evaluation_time
                .as_secs_f64(),
            event_processing_time_seconds: evaluation
                .evaluation_metadata
                .event_processing_time
                .as_secs_f64(),
        },
    })
}

fn add_contribution(
    contribution_sums: &mut BTreeMap<ContributionKey, Complex<F<f64>>>,
    key: ContributionKey,
    value: Complex<F<f64>>,
) {
    *contribution_sums.entry(key).or_insert_with(zero_complex) += value;
}

fn contribution_record(key: ContributionKey, weight: Complex<F<f64>>) -> ContributionRecord {
    ContributionRecord {
        label: contribution_label(&key),
        contribution: key.contribution,
        graph_id: key.graph_id,
        graph_name: key.graph_name,
        graph_group_id: key.graph_group_id,
        orientation_id: key.orientation_id,
        cut_id: key.cut_id,
        cut_edges: key.cut_edges,
        lmb_sample_id: key.lmb_sample_id,
        weight: complex_json(weight),
    }
}

fn contribution_label(key: &ContributionKey) -> String {
    let graph = key
        .graph_name
        .as_deref()
        .map(ToString::to_string)
        .unwrap_or_else(|| format!("#{}", key.graph_id));
    let edge_label = if key.cut_edges.is_empty() {
        "[]".to_string()
    } else {
        format!(
            "[{}]",
            key.cut_edges
                .iter()
                .map(ToString::to_string)
                .collect::<Vec<_>>()
                .join(",")
        )
    };
    let mut parts = vec![
        key.contribution.clone(),
        format!("graph={graph}"),
        format!("cut={} edges={edge_label}", key.cut_id),
    ];
    if let Some(orientation_id) = key.orientation_id {
        parts.push(format!("orientation={orientation_id}"));
    }
    if let Some(lmb_sample_id) = key.lmb_sample_id {
        parts.push(format!("lmb_sample={lmb_sample_id}"));
    }
    parts.join(" ")
}

fn additional_weight_key_label(key: AdditionalWeightKey) -> Option<String> {
    match key {
        AdditionalWeightKey::FullMultiplicativeFactor => None,
        AdditionalWeightKey::Original => Some("original".to_string()),
        AdditionalWeightKey::ThresholdCounterterm { subset_index } => {
            Some(format!("threshold_counterterm_{subset_index}"))
        }
        AdditionalWeightKey::AmplitudeThresholdCounterterm {
            esurface_id,
            overlap_group,
        } => Some(format!("ct_{esurface_id}_{overlap_group}")),
        AdditionalWeightKey::AmplitudeThresholdCountertermVariant {
            variant_id,
            esurface_id,
            overlap_group,
        } => Some(format!(
            "ct_variant_{variant_id}_{esurface_id}_{overlap_group}"
        )),
    }
}

fn zero_complex() -> Complex<F<f64>> {
    Complex::new(F(0.0), F(0.0))
}

fn complex_json(value: Complex<F<f64>>) -> ComplexJson {
    ComplexJson {
        re: value.re.0,
        im: value.im.0,
    }
}

fn approach_progress_bar(len: u64) -> Arc<ProgressBar> {
    let bar = if io::stderr().is_terminal() {
        ProgressBar::new(len)
    } else {
        ProgressBar::hidden()
    };
    bar.set_style(gammalooprs::utils::long_running_progress_style());
    Arc::new(bar)
}

fn progress_message(
    axis_index: usize,
    axis_count: usize,
    axis_point_index: usize,
    axis_point_count: usize,
    completed_total: usize,
    total: usize,
) -> String {
    format!(
        "{} {} {} {} {} {}",
        "axis".bright_cyan().bold(),
        format!("{:>2}/{:<2}", axis_index + 1, axis_count)
            .bright_white()
            .bold(),
        "point".bright_cyan().bold(),
        format!("{:>5}/{:<5}", axis_point_index + 1, axis_point_count)
            .bright_white()
            .bold(),
        "total".bright_cyan().bold(),
        format!("{:>5}/{:<5}", completed_total, total)
            .bright_green()
            .bold()
    )
}

fn relative_path_display(path: &Path) -> String {
    if path.is_absolute() {
        if let Ok(current_dir) = env::current_dir() {
            if let Ok(relative) = path.strip_prefix(current_dir) {
                return relative.display().to_string();
            }
        }
    }
    path.display().to_string()
}

fn shell_quote(value: &str) -> String {
    if value
        .chars()
        .all(|ch| ch.is_ascii_alphanumeric() || matches!(ch, '/' | '.' | '_' | '-' | '='))
    {
        value.to_string()
    } else {
        format!("'{}'", value.replace('\'', "'\\''"))
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    fn command_with_spacing(n_points: usize, logarithmic: bool, min_abs_t: f64) -> Approach {
        Approach {
            point: vec![0.5, 0.5],
            approach_axes: vec!["0.1,0.0".to_string()],
            n_points,
            logarithmic,
            min_abs_t,
            ..Default::default()
        }
    }

    #[test]
    fn linear_t_values_are_ordered_around_zero() {
        let command = command_with_spacing(2, false, 1.0e-6);
        assert_eq!(
            command.t_values(ApproachSpacing::Linear).unwrap(),
            vec![-1.0, -0.5, 0.0, 0.5, 1.0]
        );
    }

    #[test]
    fn linear_t_values_use_physical_maximum() {
        let command = Approach {
            max_abs_t: 0.2,
            ..command_with_spacing(2, false, 1.0e-6)
        };
        assert_eq!(
            command.t_values(ApproachSpacing::Linear).unwrap(),
            vec![-0.2, -0.1, 0.0, 0.1, 0.2]
        );
    }

    #[test]
    fn logarithmic_t_values_use_min_abs_t() {
        let command = command_with_spacing(3, true, 1.0e-4);
        let t_values = command.t_values(ApproachSpacing::Logarithmic).unwrap();
        assert_eq!(t_values.len(), 7);
        assert!((t_values[0] + 1.0).abs() < 1.0e-14);
        assert!((t_values[2] + 1.0e-4).abs() < 1.0e-14);
        assert_eq!(t_values[3], 0.0);
        assert!((t_values[4] - 1.0e-4).abs() < 1.0e-14);
        assert!((t_values[6] - 1.0).abs() < 1.0e-14);
    }

    #[test]
    fn logarithmic_t_values_span_physical_bounds() {
        let command = Approach {
            max_abs_t: 1.0e-2,
            ..command_with_spacing(3, true, 1.0e-6)
        };
        let t_values = command.t_values(ApproachSpacing::Logarithmic).unwrap();
        assert!((t_values[0] + 1.0e-2).abs() < 1.0e-14);
        assert!((t_values[1] + 1.0e-4).abs() < 1.0e-14);
        assert!((t_values[2] + 1.0e-6).abs() < 1.0e-14);
        assert_eq!(t_values[3], 0.0);
        assert!((t_values[4] - 1.0e-6).abs() < 1.0e-14);
        assert!((t_values[5] - 1.0e-4).abs() < 1.0e-14);
        assert!((t_values[6] - 1.0e-2).abs() < 1.0e-14);
    }

    #[test]
    fn skip_midpoint_omits_zero_from_t_values() {
        let command = Approach {
            skip_midpoint: true,
            ..command_with_spacing(2, false, 1.0e-6)
        };
        assert_eq!(
            command.t_values(ApproachSpacing::Linear).unwrap(),
            vec![-1.0, -0.5, 0.5, 1.0]
        );
    }

    #[test]
    fn missing_skip_midpoint_deserializes_to_false() {
        let mut serialized = serde_json::to_value(command_with_spacing(2, false, 1.0e-6)).unwrap();
        serialized.as_object_mut().unwrap().remove("skip_midpoint");
        let command: Approach = serde_json::from_value(serialized).unwrap();
        assert!(!command.skip_midpoint);
    }

    #[test]
    fn new_path_metadata_uses_coherent_serde_defaults() {
        let mut serialized = serde_json::to_value(command_with_spacing(2, false, 1.0e-6)).unwrap();
        let object = serialized.as_object_mut().unwrap();
        object.remove("axis_labels");
        object.remove("max_abs_t");
        object.remove("parameter_name");
        let command: Approach = serde_json::from_value(serialized).unwrap();
        assert!(command.axis_labels.is_empty());
        assert_eq!(command.max_abs_t, DEFAULT_MAX_ABS_T);
        assert_eq!(command.parameter_name, "t");
    }

    #[test]
    fn parse_axis_rejects_dimension_mismatch() {
        let err = parse_axis("1.0,0.0,2.0", 2).unwrap_err();
        assert!(format!("{err:#}").contains("Approach axis has 3 components"));
    }

    #[test]
    fn parse_axis_rejects_non_finite_components() {
        let err = parse_axis("NaN,0.0", 2).unwrap_err();
        assert!(format!("{err:#}").contains("component 0 must be finite"));
    }

    #[test]
    fn axes_receive_trimmed_user_labels() {
        let command = Approach {
            point: vec![0.5, 0.5],
            approach_axes: vec!["0.1,0.0".to_string(), "0.0,0.2".to_string()],
            axis_labels: vec![" soft ".to_string(), "threshold".to_string()],
            n_points: 1,
            ..Default::default()
        };
        command.validate_cli().unwrap();
        let axes = command.parse_axes().unwrap();
        assert_eq!(axes[0].label, "soft");
        assert_eq!(axes[0].vector, vec![0.1, 0.0]);
        assert_eq!(axes[1].label, "threshold");
        assert_eq!(axes[1].vector, vec![0.0, 0.2]);
    }

    #[test]
    fn axis_labels_must_be_complete_nonblank_and_unique() {
        for labels in [
            vec!["one".to_string()],
            vec!["one".to_string(), " ".to_string()],
            vec!["one".to_string(), "one".to_string()],
        ] {
            let command = Approach {
                point: vec![0.5, 0.5],
                approach_axes: vec!["0.1,0.0".to_string(), "0.0,0.2".to_string()],
                axis_labels: labels,
                n_points: 1,
                ..Default::default()
            };
            assert!(command.validate_cli().is_err());
        }
    }

    #[test]
    fn coordinate_jobs_skip_points_outside_unit_hypercube() {
        let command = Approach {
            point: vec![0.95, 0.5],
            approach_axes: vec!["0.1,0.0".to_string()],
            n_points: 1,
            ..Default::default()
        };
        let axes = command.parse_axes().unwrap();
        let t_values = command.t_values(ApproachSpacing::Linear).unwrap();
        let jobs = command.build_jobs(&axes, &t_values);
        assert!(jobs[2].skip_reason.is_some());
        assert!(jobs[1].skip_reason.is_none());
    }

    #[test]
    fn zero_n_cores_is_rejected() {
        let command = Approach {
            point: vec![0.5, 0.5],
            approach_axes: vec!["0.1,0.0".to_string()],
            n_points: 1,
            n_cores: Some(0),
            ..Default::default()
        };
        let err = command.validate_cli().unwrap_err();
        assert!(format!("{err:#}").contains("--n-cores must be at least 1"));
    }

    #[test]
    fn invalid_log_min_abs_t_is_rejected() {
        let command = Approach {
            point: vec![0.5, 0.5],
            approach_axes: vec!["0.1,0.0".to_string()],
            n_points: 1,
            logarithmic: true,
            min_abs_t: 0.0,
            ..Default::default()
        };
        let err = command.validate_cli().unwrap_err();
        assert!(format!("{err:#}").contains("--min-abs-t must be finite and greater than 0"));
    }

    #[test]
    fn invalid_physical_parameter_bounds_are_rejected() {
        for (min_abs_t, max_abs_t, logarithmic) in [(1.0e-2, 1.0e-3, true), (1.0e-3, 0.0, false)] {
            let command = Approach {
                point: vec![0.5, 0.5],
                approach_axes: vec!["0.1,0.0".to_string()],
                n_points: 1,
                min_abs_t,
                max_abs_t,
                logarithmic,
                ..Default::default()
            };
            assert!(command.validate_cli().is_err());
        }
    }

    #[test]
    fn non_finite_midpoint_and_blank_parameter_name_are_rejected() {
        let non_finite = Approach {
            point: vec![f64::INFINITY, 0.5],
            approach_axes: vec!["0.1,0.0".to_string()],
            n_points: 1,
            ..Default::default()
        };
        assert!(format!("{:#}", non_finite.validate_cli().unwrap_err()).contains("must be finite"));

        let blank_parameter = Approach {
            point: vec![0.5, 0.5],
            approach_axes: vec!["0.1,0.0".to_string()],
            n_points: 1,
            parameter_name: " ".to_string(),
            ..Default::default()
        };
        assert!(format!("{:#}", blank_parameter.validate_cli().unwrap_err())
            .contains("--parameter-name must not be blank"));
    }

    #[test]
    fn threshold_summary_is_scoped_and_grouped_by_graph_component() {
        use gammalooprs::observables::events::{
            GenericThresholdCountertermComponentWeight, ThresholdCountertermComponentOccurrence,
        };

        let component = |component_id, weighted, evaluation_skipped| {
            GenericThresholdCountertermComponentWeight {
                component_id,
                occurrence: ThresholdCountertermComponentOccurrence::Amplitude {
                    raised_esurface_id: component_id,
                    overlap_group: 0,
                },
                multiplier_values: Default::default(),
                effective_multiplier: F(1.0),
                bare: None,
                weighted: Complex::new_re(F(weighted)),
                evaluation_skipped,
            }
        };
        let mut accumulator = ThresholdCountertermSummaryAccumulator::default();
        let first = GenericThresholdCountertermEventInfo {
            original: Complex::new_re(F(10.0)),
            components: vec![component(2, 1.0, false), component(2, 2.0, true)],
        };
        let second = GenericThresholdCountertermEventInfo {
            original: Complex::new_re(F(20.0)),
            components: vec![component(3, -4.0, false)],
        };
        let third = GenericThresholdCountertermEventInfo {
            original: Complex::new_re(F(5.0)),
            components: vec![component(2, 1.0, false)],
        };
        accumulator.add(1, &Complex::new_re(F(13.0)), &first);
        accumulator.add(1, &Complex::new_re(F(16.0)), &second);
        accumulator.add(2, &Complex::new_re(F(6.0)), &third);

        let summary = accumulator.into_json().unwrap();
        assert_eq!(summary.original.re, 35.0);
        assert_eq!(summary.counterterm_sum.re, 0.0);
        assert_eq!(summary.reconstructed_total.re, 35.0);
        assert_eq!(summary.decomposed_event_sum.re, 35.0);
        assert_eq!(summary.closure.re, 0.0);
        assert_eq!(summary.components.len(), 3);
        assert_eq!(
            (
                summary.components[0].graph_id,
                summary.components[0].component_id,
                summary.components[0].weighted.re,
                summary.components[0].occurrence_count,
                summary.components[0].skipped_count,
            ),
            (1, 2, 3.0, 2, 1),
        );
        assert_eq!(
            (
                summary.components[2].graph_id,
                summary.components[2].component_id,
            ),
            (2, 2),
        );
    }

    #[test]
    fn multiplicative_factor_is_not_an_additive_weight_label() {
        assert_eq!(
            additional_weight_key_label(AdditionalWeightKey::FullMultiplicativeFactor),
            None,
        );
        assert_eq!(
            additional_weight_key_label(AdditionalWeightKey::Original).as_deref(),
            Some("original"),
        );
    }

    #[test]
    fn simulated_parallel_completion_is_sorted_back_to_job_order() {
        let command = Approach {
            point: vec![0.5, 0.5],
            approach_axes: vec!["0.1,0.0".to_string(), "0.0,0.2".to_string()],
            n_points: 2,
            ..Default::default()
        };
        let axes = command.parse_axes().unwrap();
        let t_values = command.t_values(ApproachSpacing::Linear).unwrap();
        let jobs = command.build_jobs(&axes, &t_values);
        let mut records = (0..2)
            .flat_map(|worker_index| jobs.iter().skip(worker_index).step_by(2))
            .map(|job| ApproachPointRecord {
                index: job.index,
                axis_index: job.axis_index,
                axis_point_index: job.axis_point_index,
                t: job.t,
                point: job.point.clone(),
                status: "skipped",
                skip_reason: None,
                evaluation: None,
            })
            .collect::<Vec<_>>();

        assert_ne!(
            records
                .iter()
                .map(|record| record.index)
                .collect::<Vec<_>>(),
            (0..records.len()).collect::<Vec<_>>(),
            "strided worker completion should be out of original order in this simulation",
        );
        records.sort_by_key(|record| record.index);
        assert_eq!(
            records
                .iter()
                .map(|record| (record.index, record.axis_index, record.axis_point_index))
                .collect::<Vec<_>>(),
            jobs.iter()
                .map(|job| (job.index, job.axis_index, job.axis_point_index))
                .collect::<Vec<_>>(),
        );
    }
}