From d6d33b1c73b32ca952c08f890bbae3e4e58410b9 Mon Sep 17 00:00:00 2001 From: ApsisRS contributors Date: Wed, 30 Sep 2026 17:53:48 +0800 Subject: [PATCH] bench: add a parallel scale benchmark Splits the three heavy paths across scoped threads and reports the speedup and the efficiency against the thread count, with the serial time taken as the fastest of several runs so a scheduling hiccup cannot inflate the reference every number is divided by. Constellation propagation and access search are embarrassingly parallel and scale as expected until the thread count passes eight, where the machine stops gaining. Topology construction does not scale at all, and the split line shows why: the geometric visibility test is two milliseconds while folding the visible pairs into the graph is four hundred, because each insertion scans the edges already stored. --- .github/workflows/ci.yml | 1 + crates/apsis-rs/Cargo.toml | 4 + crates/apsis-rs/examples/bench_parallel.rs | 421 +++++++++++++++++++++ 3 files changed, 426 insertions(+) create mode 100644 crates/apsis-rs/examples/bench_parallel.rs diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 250d86d..0c2b240 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -41,6 +41,7 @@ jobs: run: | cargo build -p apsis-rs --features full --example bench cargo build -p apsis-rs --features full --example bench_scale + cargo build -p apsis-rs --features full --example bench_parallel - name: Documentation # Catches a doc link that no longer resolves, which is how a doc # comment drifts away from the code it describes. diff --git a/crates/apsis-rs/Cargo.toml b/crates/apsis-rs/Cargo.toml index 6b67ea9..296b695 100644 --- a/crates/apsis-rs/Cargo.toml +++ b/crates/apsis-rs/Cargo.toml @@ -221,3 +221,7 @@ required-features = ["orbit"] [[example]] name = "bench_scale" required-features = ["full"] + +[[example]] +name = "bench_parallel" +required-features = ["full"] diff --git a/crates/apsis-rs/examples/bench_parallel.rs b/crates/apsis-rs/examples/bench_parallel.rs new file mode 100644 index 0000000..bb790b5 --- /dev/null +++ b/crates/apsis-rs/examples/bench_parallel.rs @@ -0,0 +1,421 @@ +//! Parallel scale benchmark: how well the hot paths use the cores. +//! +//! The serial `bench_scale` example measures the cost of a path. This one +//! asks the next question: when the work is spread over the cores, how much +//! of the wall clock does that save, and where does the parallel version +//! stop being worth it. +//! +//! The paths chosen are the ones a real run parallelises, because each is +//! either trivially parallel or has a clear serial tail: +//! +//! - **constellation propagation** propagates every satellite independently, +//! so it is embarrassingly parallel and should scale nearly linearly. +//! - **access search** searches each track independently, likewise. +//! - **topology construction** decides each pair of satellites independently, +//! but folding the results into one graph is serial, so it is the case +//! that shows the cost of a serial tail. +//! +//! The timing compares the same total work done serially and across a range +//! of thread counts, and reports the speedup against both the serial time +//! and the ideal linear time. The two differ, and the difference is the +//! parallel overhead: thread spawn, work imbalance, and the serial tail. +//! +//! There is no warm-up, for the reason given in `bench_scale`: the first +//! call pays for whatever the path allocates on first use, and a caller +//! waits for that too. The serial time and the parallel time are both +//! measured cold, so the comparison is between two numbers the caller could +//! actually observe. +//! +//! Usage: +//! +//! ```text +//! cargo run --release -p apsis-rs --features full --example bench_parallel +//! ``` + +use apsis_access::{AccessFinder, OrbitSample}; +use apsis_coord::Eci; +use apsis_coverage::GroundPoint; +use apsis_frames::EarthRotationAngle; +use apsis_geodesy::Geodetic; +use apsis_isl::mutual_visibility; +use apsis_math::Vec3; +use apsis_sgp4::{GravityModel, Satellite}; +use apsis_tle::TwoLineElement; +use apsis_topology::Topology; +use std::hint::black_box; +use std::time::{Duration, Instant}; + +/// The element set used as the source of every propagated track. +const ISS: &str = concat!( + "1 25544U 98067A 08264.51782528 -.00002182 00000-0 -11606-4 0 2927\r\n", + "2 25544 51.6416 247.4627 0006703 130.5360 325.0288 15.72125391563537\r\n", +); + +/// The body radius, in kilometers, taken from the shared constants. +const BODY_RADIUS_KM: f64 = apsis_core::constant::EARTH_EQUATORIAL_RADIUS; + +/// The sample step, in seconds. +const TRACK_STEP_SECONDS: f64 = 60.0; + +/// Runs a path once and returns how long it took. +fn measure(path: F) -> Duration +where + F: FnOnce() -> T, +{ + let start = Instant::now(); + black_box(path()); + start.elapsed() +} + +/// Runs a path `repeats` times and returns the fastest duration. +/// +/// The fastest of several runs is the least contaminated by whatever else +/// the machine is doing, and it is the figure that scales cleanly, because a +/// scheduling hiccup on one run cannot inflate the reference every later +/// number is divided by. The cold cost of the very first call is still +/// reported on its own, so nothing is hidden. +fn measure_best(repeats: u32, path: F) -> Duration +where + F: Fn() -> T, +{ + let mut best = Duration::MAX; + for _ in 0..repeats { + let elapsed = measure(&path); + if elapsed < best { + best = elapsed; + } + } + best +} + +/// Prints one row. +/// +/// The efficiency is the speedup divided by the thread count, so one means +/// the path used every core fully. A value well below one means the threads +/// were idle, which is what a serial tail looks like. +fn row(threads: usize, wall: Duration, work: Duration, cores: usize) { + let wall_ms = wall.as_secs_f64() * 1e3; + let work_ms = work.as_secs_f64() * 1e3; + let speedup = work_ms / wall_ms; + let efficiency = speedup / threads as f64; + println!( + "{threads:>6} {wall_ms:>11.3} {work_ms:>11.3} {speedup:>9.2} {efficiency:>9.2} {cores:>8}" + ); +} + +/// Returns a prepared satellite. +fn satellite() -> Satellite { + let element = TwoLineElement::parse(ISS).expect("the sample parses"); + Satellite::from_element(&element, GravityModel::WGS72).expect("the orbit propagates") +} + +/// The epoch of the sample element set, as a Modified Julian Date. +fn epoch_mjd() -> f64 { + let element = TwoLineElement::parse(ISS).expect("the sample parses"); + element + .epoch + .mjd(apsis_time::TimeScale::Utc) + .expect("the epoch is inside the leap second table") +} + +/// Returns a track of Earth-fixed samples for one satellite. +fn track(satellite: &Satellite, samples: usize, epoch_mjd: f64) -> Vec { + (0..samples) + .map(|index| { + let seconds = index as f64 * TRACK_STEP_SECONDS; + let state = satellite + .propagate_minutes(seconds / 60.0) + .expect("the track propagates"); + let angle = + EarthRotationAngle::from_ut1_days(epoch_mjd - 51_544.5 + seconds / 86_400.0); + let fixed = Eci::new(state.position).to_ecef(angle).vector(); + OrbitSample::new(seconds, fixed) + }) + .collect() +} + +/// Propagates a whole constellation over a day, serially. +/// +/// The serial time is the reference every parallel number is divided by. +fn propagate_serial(satellite: &Satellite, satellites: usize, samples: usize) -> usize { + let mut count = 0_usize; + for _ in 0..satellites { + for index in 0..samples { + let seconds = index as f64 * TRACK_STEP_SECONDS; + let state = satellite + .propagate_minutes(seconds / 60.0) + .expect("the track propagates"); + black_box(state.position); + count += 1; + } + } + count +} + +/// Propagates a whole constellation over a day across `threads` threads. +/// +/// Each thread takes a disjoint slice of the satellites, so the work splits +/// evenly and no two threads write the same state. The `scope` form of +/// spawning borrows the satellite rather than requiring it to outlive the +/// threads, which is what lets the reference be shared without cloning. +fn propagate_parallel( + satellite: &Satellite, + satellites: usize, + samples: usize, + threads: usize, +) -> usize { + let chunk = satellites.div_ceil(threads); + let total = std::sync::atomic::AtomicUsize::new(0); + std::thread::scope(|scope| { + for worker in 0..threads { + let first = worker * chunk; + if first >= satellites { + break; + } + let last = (first + chunk).min(satellites); + let total = &total; + scope.spawn(move || { + let mut count = 0_usize; + for _ in first..last { + for index in 0..samples { + let seconds = index as f64 * TRACK_STEP_SECONDS; + let state = satellite + .propagate_minutes(seconds / 60.0) + .expect("the track propagates"); + black_box(state.position); + count += 1; + } + } + total.fetch_add(count, std::sync::atomic::Ordering::Relaxed); + }); + } + }); + total.load(std::sync::atomic::Ordering::Relaxed) +} + +/// Searches access for every track, serially. +fn access_serial(tracks: &[Vec], finder: &AccessFinder) -> usize { + let mut passes = 0_usize; + for track in tracks { + passes += finder.find(track).expect("the track is valid").count(); + } + passes +} + +/// Searches access for every track across `threads` threads. +/// +/// The tracks are read only and the finder is a copy, so each thread can +/// take a disjoint slice without synchronisation. +fn access_parallel(tracks: &[Vec], finder: &AccessFinder, threads: usize) -> usize { + let chunk = tracks.len().div_ceil(threads); + let total = std::sync::atomic::AtomicUsize::new(0); + std::thread::scope(|scope| { + for worker in 0..threads { + let first = worker * chunk; + if first >= tracks.len() { + break; + } + let last = (first + chunk).min(tracks.len()); + let total = &total; + let finder = *finder; + scope.spawn(move || { + let mut passes = 0_usize; + for track in &tracks[first..last] { + passes += finder.find(track).expect("the track is valid").count(); + } + total.fetch_add(passes, std::sync::atomic::Ordering::Relaxed); + }); + } + }); + total.load(std::sync::atomic::Ordering::Relaxed) +} + +/// Builds a topology serially, deciding each pair and folding it in. +fn topology_serial(positions: &[Vec3], satellites: usize, range: f64) -> usize { + let mut topology = Topology::new(satellites); + for first in 0..satellites { + for second in 0..satellites { + if first != second + && mutual_visibility(positions[first], positions[second], range) + .is_some_and(|visible| visible) + { + let _ = topology.add_edge(first, second, 1.0); + } + } + } + topology.edge_count() +} + +/// Builds a topology across `threads` threads. +/// +/// The geometry of each pair is independent, so each thread computes the +/// visible pairs in its own slice and returns them. Folding them into the +/// single graph is serial, which is the point of measuring this path: it +/// shows how much of the time the parallel part actually saves. +fn topology_parallel(positions: &[Vec3], satellites: usize, range: f64, threads: usize) -> usize { + let chunk = satellites.div_ceil(threads); + let results: Vec> = std::thread::scope(|scope| { + let handles: Vec<_> = (0..threads) + .map(|worker| { + let first = worker * chunk; + let last = (first + chunk).min(satellites); + scope.spawn(move || { + let mut pairs = Vec::new(); + for first in first..last { + for second in 0..satellites { + if first != second + && mutual_visibility(positions[first], positions[second], range) + .is_some_and(|visible| visible) + { + pairs.push((first, second)); + } + } + } + pairs + }) + }) + .collect(); + handles + .into_iter() + .map(|handle| handle.join().expect("the worker did not panic")) + .collect() + }); + + // The serial tail: fold every visible pair into one graph. + let mut topology = Topology::new(satellites); + for pairs in &results { + for &(first, second) in pairs { + let _ = topology.add_edge(first, second, 1.0); + } + } + topology.edge_count() +} + +/// Splits topology construction into its geometric and its folding phase. +/// +/// Returns the time of the visibility test and the time of folding the +/// visible pairs into a graph, both in milliseconds. The first is +/// embarrassingly parallel, the second is a serial pass whose cost grows +/// with the square of the edge count because each insertion checks for a +/// duplicate by scanning the edges already stored. +fn topology_phase_split(positions: &[Vec3], satellites: usize, range: f64) -> (f64, f64) { + let start = Instant::now(); + let mut pairs = Vec::new(); + for first in 0..satellites { + for second in 0..satellites { + if first != second + && mutual_visibility(positions[first], positions[second], range) + .is_some_and(|visible| visible) + { + pairs.push((first, second)); + } + } + } + let visible_ms = start.elapsed().as_secs_f64() * 1e3; + + let start = Instant::now(); + let mut topology = Topology::new(satellites); + for &(first, second) in &pairs { + let _ = topology.add_edge(first, second, 1.0); + } + let fold_ms = start.elapsed().as_secs_f64() * 1e3; + + (visible_ms, fold_ms) +} + +fn main() { + let cores = std::thread::available_parallelism() + .map(|count| count.get()) + .unwrap_or(1); + println!("ApsisRS parallel scale benchmark: release build, no external dependency"); + println!("available parallelism: {cores} core(s)"); + println!("no warm-up: the cold cost of the first call is reported, not discarded"); + + let satellite = satellite(); + let epoch = epoch_mjd(); + + // The thread counts swept: one is the serial reference, then the powers + // of two up to the machine, then the machine's own count. + let mut thread_counts = vec![1_usize, 2, 4, 8]; + if cores > 8 { + thread_counts.push(cores); + } + thread_counts.retain(|count| *count <= cores.max(2)); + thread_counts.dedup(); + + // ---- One: constellation propagation over a day. ---------------------- + { + let samples = 1_441_usize; // one day at a one minute step + for &satellites in &[550_usize, 1_000, 2_000] { + println!( + "\none: propagate {satellites} satellite(s) over 1 day at 60 s ({samples} samples each)" + ); + let serial = measure_best(3, || propagate_serial(&satellite, satellites, samples)); + let work = serial; + for &threads in &thread_counts { + let wall = measure_best(3, || { + propagate_parallel(&satellite, satellites, samples, threads) + }); + row(threads, wall, work, cores); + } + } + } + + // ---- Two: access search. --------------------------------------------- + { + for &satellites in &[64_usize, 256] { + let tracks: Vec> = (0..satellites) + .map(|_| track(&satellite, 120, epoch)) + .collect(); + let finder = AccessFinder { + station: GroundPoint::new(Geodetic::new(0.0, 0.0, 0.0), 10.0_f64.to_radians()), + sample_step_seconds: TRACK_STEP_SECONDS, + crossing_tolerance_seconds: 0.01, + }; + println!("\ntwo: access search over {satellites} track(s), 2 h each"); + let serial = measure_best(5, || access_serial(&tracks, &finder)); + let work = serial; + for &threads in &thread_counts { + let wall = measure_best(5, || access_parallel(&tracks, &finder, threads)); + row(threads, wall, work, cores); + } + } + } + + // ---- Three: topology construction. ------------------------------------ + { + for &satellites in &[120_usize, 300] { + let positions: Vec = (0..satellites) + .map(|index| { + let angle = index as f64 * 0.4; + Vec3::new( + (BODY_RADIUS_KM + 550.0) * angle.cos(), + (BODY_RADIUS_KM + 550.0) * angle.sin(), + index as f64 * 20.0, + ) + }) + .collect(); + println!("\nthree: build topology over {satellites} satellite(s), 4500 km range"); + let serial = measure_best(2, || topology_serial(&positions, satellites, 4_500.0)); + // The split shows where the time actually goes: the geometric + // visibility test is parallel, but folding the visible pairs into + // the graph is a serial pass whose cost is superlinear in the + // edge count. Reporting the two separately is the difference + // between "parallelism did not help" and knowing why. + let (vis_ms, fold_ms) = topology_phase_split(&positions, satellites, 4_500.0); + println!( + " split: visibility {vis_ms:.3} ms (parallel), edge fold {fold_ms:.3} ms (serial tail)" + ); + for &threads in &thread_counts { + let wall = measure_best(2, || { + topology_parallel(&positions, satellites, 4_500.0, threads) + }); + row(threads, wall, serial, cores); + } + } + } + + println!(); + println!("speedup is serial time over wall time; efficiency is speedup over thread count"); + println!("an efficiency near one means the cores were kept busy; a low value is a serial tail"); +}