diff --git a/CHANGELOG.md b/CHANGELOG.md index 830a68b8..f725920c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,14 @@ removed no sooner than the next major (see `docs/API_STABILITY.md`). ### Changed +- Improve resident trimmed ICP moment precision with a default DoubleSingle + route: integer product residuals and preserved f32 rounding carry two words + through coordinate normalization, centroid division and covariance products. + Explicit `GpuIcpTrimPrecision::CompensatedF32` retains the earlier arithmetic. + Fraction one keeps canonical untrimmed results; matching, all-gated fitness, + CPU Horn solving and fixed 176-byte readbacks remain. Large gated registration + regressions now enforce the unchanged CPU pose-agreement tolerance. + - Add explicit resident GPU ICP trimming under `register-icp-gpu`. Stable distance/source-index selection retains floor(trim_fraction*gated_count) pairs for moments while all gated pairs still define fitness/stopping. diff --git a/crates/spatialrust-gpu/src/kernels/icp_reduction.rs b/crates/spatialrust-gpu/src/kernels/icp_reduction.rs index fb6cfb59..eeb868d6 100644 --- a/crates/spatialrust-gpu/src/kernels/icp_reduction.rs +++ b/crates/spatialrust-gpu/src/kernels/icp_reduction.rs @@ -165,12 +165,19 @@ impl GpuIcpInput { transform: Isometry3, max_distance: f32, ) -> SpatialResult { - self.moments_impl(runtime, transform, max_distance, None) + self.moments_impl( + runtime, + transform, + max_distance, + None, + trim::GpuIcpTrimPrecision::DoubleSingle, + ) } /// Submits stable GPU selection of the lowest-distance floor(fraction*N) /// gated pairs, breaking distance ties by original source index. Moments - /// use retained pairs; fitness uses ALL gated pairs. Binary64 fraction + /// use retained pairs with DoubleSingle normalized arithmetic; fitness + /// uses ALL gated pairs. Binary64 fraction /// multiplication/floor is emulated with integer limbs on the device. /// No match/point payload readback or host wait occurs. A fraction of one /// delegates to the canonical untrimmed kernels; no CPU fallback occurs. @@ -180,6 +187,28 @@ impl GpuIcpInput { transform: Isometry3, max_distance: f32, fraction: f64, + ) -> SpatialResult { + self.trimmed_moments_resident_with_precision( + runtime, + transform, + max_distance, + fraction, + trim::GpuIcpTrimPrecision::DoubleSingle, + ) + } + + /// Explicit choice of retained-moment arithmetic. DoubleSingle carries + /// two f32 words through normalized offsets, division and products; f32 + /// transforms/search stay unchanged. CompensatedF32 retains the earlier + /// coordinate/product precision. Fraction one uses canonical ordinary + /// moments for either choice. Transfers/ownership/status rules are identical. + pub fn trimmed_moments_resident_with_precision( + &self, + runtime: &WgpuRuntime, + transform: Isometry3, + max_distance: f32, + fraction: f64, + precision: trim::GpuIcpTrimPrecision, ) -> SpatialResult { if !fraction.is_finite() || fraction <= 0.0 || fraction > 1.0 { return Err(SpatialError::InvalidArgument( @@ -187,8 +216,13 @@ impl GpuIcpInput { )); } let trimmed = fraction < 1.0; - let inner = - self.moments_impl(runtime, transform, max_distance, trimmed.then_some(fraction))?; + let inner = self.moments_impl( + runtime, + transform, + max_distance, + trimmed.then_some(fraction), + precision, + )?; Ok(GpuIcpTrimmedMoments { inner, trimmed }) } @@ -198,6 +232,7 @@ impl GpuIcpInput { transform: Isometry3, max_distance: f32, fraction: Option, + precision: trim::GpuIcpTrimPrecision, ) -> SpatialResult { let q = transform.rotation(); let t = transform.translation(); @@ -266,10 +301,10 @@ impl GpuIcpInput { let partials = scratch.acquire("icp-partials", u64::from(groups) * 144); let selection = fraction.map(|_| scratch.acquire("icp-trim-selection", 32)); let cache = if fraction.is_some() { - &runtime + runtime .icp_trim_pipelines .get_or_init(|| trim::TrimPipelines::new(runtime.device())) - .reduction + .reduction(runtime.device(), precision) } else { runtime .icp_reduction_pipelines diff --git a/crates/spatialrust-gpu/src/kernels/icp_trim.rs b/crates/spatialrust-gpu/src/kernels/icp_trim.rs index de1bae62..fdc4ba9a 100644 --- a/crates/spatialrust-gpu/src/kernels/icp_trim.rs +++ b/crates/spatialrust-gpu/src/kernels/icp_trim.rs @@ -3,6 +3,19 @@ use super::*; pub(super) const SHADER: &str = include_str!("../shaders/icp_trim.wgsl"); +/// Arithmetic used for retained ICP moments; transform/search remain f32. +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +pub enum GpuIcpTrimPrecision { + /// Canonical f32 normalized coordinates/products with compensated sums. + /// Lower arithmetic cost, but coordinate/product rounding can bias moments. + CompensatedF32, + /// Two f32 words through offsets, normalization, mean division and products. + /// Carries quotient/product residuals to reduce moment bias; hardware f32 + /// transforms/search and their gate/tie precision remain unchanged. + #[default] + DoubleSingle, +} + #[repr(C)] #[derive(Clone, Copy, Pod, Zeroable)] pub(super) struct TrimParams { @@ -48,7 +61,8 @@ impl TrimParams { } pub(crate) struct TrimPipelines { - pub(super) reduction: IcpReductionPipelines, + reduction: std::sync::OnceLock, + precise: std::sync::OnceLock, layout: wgpu::BindGroupLayout, histogram: wgpu::ComputePipeline, select: wgpu::ComputePipeline, @@ -96,13 +110,29 @@ impl TrimPipelines { }) }); Self { - reduction: IcpReductionPipelines::with_shader(device, &trimmed_shader(), true), + reduction: std::sync::OnceLock::new(), + precise: std::sync::OnceLock::new(), layout, histogram, select, } } + pub(super) fn reduction( + &self, + device: &wgpu::Device, + precision: GpuIcpTrimPrecision, + ) -> &IcpReductionPipelines { + match precision { + GpuIcpTrimPrecision::CompensatedF32 => self.reduction.get_or_init(|| { + IcpReductionPipelines::with_shader(device, &trimmed_shader(), true) + }), + GpuIcpTrimPrecision::DoubleSingle => self.precise.get_or_init(|| { + IcpReductionPipelines::with_shader(device, &precise_shader(), true) + }), + } + } + pub(super) fn encode( &self, runtime: &WgpuRuntime, @@ -192,6 +222,64 @@ fn trimmed_shader() -> String { shader } +fn precise_shader() -> String { + let mut shader = trimmed_shader(); + shader.push_str(include_str!("../shaders/icp_trim_moments.wgsl")); + replace_once(&mut shader, "var value=0.0;", "var value=vec2(0.0);"); + replace_once( + &mut shader, + "value=((p-anchor)/params.scale)[term];", + "value=ds_offset(p[term],anchor[term],params.scale);", + ); + replace_once( + &mut shader, + "value=((q-params.target_anchor.xyz)/params.scale)[term-3u];", + "value=ds_offset(q[term-3u],params.target_anchor[term-3u],params.scale);", + ); + replace_once(&mut shader, "value=dot3(d,d);", "value=vec2(dot3(d,d),0.0);"); + replace_once(&mut shader, "if !finite(value) {", "if !finite(value.x) || !finite(value.y) {"); + replace_once( + &mut shader, + "total=add(total,value);", + "if term==15u { total=add(total,value.x); } else { total=precise_merge(total,value); }", + ); + replace_once(&mut shader, "result.values[term]=merge(vec2(total.x/f32(count),0.0),vec2(total.y/f32(count),0.0));", "result.values[term]=ds_divide(total,f32(count));"); + replace_once( + &mut shader, + "let mean_s=result.values[row].x+result.values[row].y;", + "let mean_s=result.values[row];", + ); + replace_once( + &mut shader, + "let mean_t=result.values[col+3u].x+result.values[col+3u].y;", + "let mean_t=result.values[col+3u];", + ); + replace_once( + &mut shader, + "let s=(queries[i].xyz-result.source_anchor.xyz)/params.scale;", + "let s=ds_offset(queries[i][row],result.source_anchor[row],params.scale);", + ); + replace_once(&mut shader, "let t=(target_points[matches[i].index].xyz-result.target_anchor.xyz)/params.scale;", "let t=ds_offset(target_points[matches[i].index][col],result.target_anchor[col],params.scale);"); + replace_once( + &mut shader, + "total=add(total,(s[row]-mean_s)*(t[col]-mean_t));", + "total=precise_merge(total,ds_product(precise_merge(s,-mean_s),precise_merge(t,-mean_t)));", + ); + replace_once(&mut shader, "if lane = vec![ + [1.0, 0.0, 3.0, 0.0], + [1.0, 2.0f32.powi(-24), 1.0, -2.0f32.powi(-24)], + [16777216.0, 0.0, 1.0, 0.0], + [8.2388, 1e-7, 209642.0, 0.0], + [-2548.44, 0.00007, 209642.0, 0.0], + [-64.0, 0.0, 128.7, 0.0], + [-0.78, 1e-8, 131.25, 1e-7], + ]; + let mut seed = 20261013u32; + for i in 0..512 { + seed = seed.wrapping_mul(1664525).wrapping_add(1013904223); + let a = f32::from_bits( + ((100 + (i % 50)) as u32) << 23 | (seed & 0x7fffff) | ((i as u32 & 1) << 31), + ); + seed = seed.wrapping_mul(1664525).wrapping_add(1013904223); + let b = f32::from_bits( + ((100 + ((i * 7) % 50)) as u32) << 23 + | (seed & 0x7fffff) + | (((i as u32 >> 1) & 1) << 31), + ); + inputs.push([a, a * 2.0f32.powi(-27), b, -b * 2.0f32.powi(-27)]); + } + let source = include_str!("../shaders/icp_reduction.wgsl"); + let start = source.find("fn add(").unwrap(); + let end = source.find("fn valid_match(").unwrap(); + let shader=format!("struct Params {{ unused:u32,count:u32,pad0:u32,pad1:u32 }}\n@group(0) @binding(0) var params:Params;\n@group(0) @binding(1) var inputs:array>;\n@group(0) @binding(2) var output:array>;\n{}\n{}\n@compute @workgroup_size(64) fn test_math(@builtin(global_invocation_id) id:vec3) {{ if id.x>=params.count {{ return; }} let v=inputs[id.x]; output[4u*id.x]=precise_add(v.xy,v.z); output[4u*id.x+1u]=two_product(v.x,v.z); output[4u*id.x+2u]=ds_divide(v.xy,v.z); output[4u*id.x+3u]=ds_product(v.xy,v.zw); }}", &source[start..end],include_str!("../shaders/icp_trim_moments.wgsl")); + let module = runtime.device().create_shader_module(wgpu::ShaderModuleDescriptor { + label: None, + source: wgpu::ShaderSource::Wgsl(shader.into()), + }); + let pipeline = runtime.device().create_compute_pipeline(&wgpu::ComputePipelineDescriptor { + label: None, + layout: None, + module: &module, + entry_point: Some("test_math"), + compilation_options: Default::default(), + cache: None, + }); + let uniform = runtime.device().create_buffer_init(&wgpu::util::BufferInitDescriptor { + label: None, + contents: bytemuck::cast_slice(&[0u32, inputs.len() as u32, 0, 0]), + usage: wgpu::BufferUsages::UNIFORM, + }); + let input = runtime.device().create_buffer_init(&wgpu::util::BufferInitDescriptor { + label: None, + contents: bytemuck::cast_slice(&inputs), + usage: wgpu::BufferUsages::STORAGE, + }); + let output = runtime.device().create_buffer(&wgpu::BufferDescriptor { + label: None, + size: inputs.len() as u64 * 32, + usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC, + mapped_at_creation: false, + }); + let entries: Vec<_> = [&uniform, &input, &output] + .iter() + .enumerate() + .map(|(binding, b)| wgpu::BindGroupEntry { + binding: binding as u32, + resource: b.as_entire_binding(), + }) + .collect(); + let group = runtime.device().create_bind_group(&wgpu::BindGroupDescriptor { + label: None, + layout: &pipeline.get_bind_group_layout(0), + entries: &entries, + }); + let mut encoder = runtime.device().create_command_encoder(&Default::default()); + dispatch(&mut encoder, &pipeline, &group, (inputs.len() as u32).div_ceil(64), 1); + runtime.queue().submit(Some(encoder.finish())); + let words = read(&runtime, &output, output.size()); + let mut maximum_relative_error = 0.0_f64; + for (i, v) in inputs.iter().enumerate() { + let a = f64::from(v[0]) + f64::from(v[1]); + let b = f64::from(v[2]) + f64::from(v[3]); + let expected = + [a + f64::from(v[2]), f64::from(v[0]) * f64::from(v[2]), a / f64::from(v[2]), a * b]; + for (op, expected) in expected.into_iter().enumerate() { + let hi = f32::from_bits(words[8 * i + 2 * op]); + let lo = f32::from_bits(words[8 * i + 2 * op + 1]); + let actual = f64::from(hi) + f64::from(lo); + maximum_relative_error = + maximum_relative_error.max((actual - expected).abs() / expected.abs().max(1e-12)); + assert!( + (actual - expected).abs() <= 1e-11 * expected.abs().max(1e-12), + "input={i} operation={op}" + ); + } + } + eprintln!( + "GPU_ICP_PRECISE_MATH comparisons={} maximum_relative_error={maximum_relative_error}", + inputs.len() * 4 + ); +} diff --git a/crates/spatialrust-gpu/src/kernels/icp_trim_precision_tests.rs b/crates/spatialrust-gpu/src/kernels/icp_trim_precision_tests.rs new file mode 100644 index 00000000..b7ab3d1e --- /dev/null +++ b/crates/spatialrust-gpu/src/kernels/icp_trim_precision_tests.rs @@ -0,0 +1,255 @@ +//! Fixed-pose independent f64 oracles; copies occur before scratch recycling. +use super::*; + +fn snapshot( + runtime: &WgpuRuntime, + input: &GpuIcpInput, + pose: Isometry3, + precision: GpuIcpTrimPrecision, +) -> (GpuIcpTrimmedStatistics, Vec<[f32; 4]>, Vec<[u32; 2]>, Vec) { + let count = input.source_count(); + let groups = plan_trim(u64::from(count), &runtime.device().limits()).unwrap(); + let r = pose.rotation().to_mat3(); + let t = pose.translation(); + let params = Params { + rows: [ + [r.m[0][0], r.m[0][1], r.m[0][2], t.x], + [r.m[1][0], r.m[1][1], r.m[1][2], t.y], + [r.m[2][0], r.m[2][1], r.m[2][2], t.z], + ], + source_anchor: input.source_anchor, + target_anchor: input.target_anchor, + count, + targets: input.target.query_count(), + groups, + unused: 0, + bound: 0.0625, + max_distance: 0.25, + scale: input.scale, + padding: 0.0, + }; + let result = GpuIcpTrimmedMoments { + inner: GpuIcpMoments { + buffer: runtime.device().create_buffer(&wgpu::BufferDescriptor { + label: None, + size: 176, + usage: wgpu::BufferUsages::STORAGE | wgpu::BufferUsages::COPY_SRC, + mapped_at_creation: false, + }), + count, + scale: input.scale, + bound: 0.0625, + identity: runtime.profiling_identity(), + }, + trimmed: true, + }; + let uniform = runtime.device().create_buffer_init(&wgpu::util::BufferInitDescriptor { + label: None, + contents: bytemuck::bytes_of(¶ms), + usage: wgpu::BufferUsages::UNIFORM, + }); + let mut scratch = ScratchBuffers::new(runtime); + let queries = scratch.acquire("precision-query", u64::from(count) * 16); + let matches = scratch.acquire("precision-matches", u64::from(count) * 8); + let partials = scratch.acquire("precision-partials", u64::from(groups) * 144); + let selection = scratch.acquire("precision-selection", 32); + let cache = runtime.icp_trim_pipelines.get_or_init(|| TrimPipelines::new(runtime.device())); + let reduction = cache.reduction(runtime.device(), precision); + let entries: Vec<_> = [ + &uniform, + input.source.buffer().unwrap(), + input.target.buffer().unwrap(), + &queries, + &matches, + &partials, + result.buffer(), + &selection, + ] + .iter() + .enumerate() + .map(|(binding, buffer)| wgpu::BindGroupEntry { + binding: binding as u32, + resource: buffer.as_entire_binding(), + }) + .collect(); + let group = runtime.device().create_bind_group(&wgpu::BindGroupDescriptor { + label: None, + layout: &reduction.layout, + entries: &entries, + }); + let mut encoder = runtime.device().create_command_encoder(&Default::default()); + dispatch(&mut encoder, &reduction.pipelines[0], &group, count.div_ceil(64), 1); + input.tree.encode_nearest(runtime, &mut encoder, &queries, count, 0.0625, &matches).unwrap(); + cache.encode( + runtime, + &mut encoder, + &mut scratch, + TrimParams::new(count, groups, 0.8, 0.0625, params.targets), + &matches, + &selection, + ); + for (index, x, y) in [(1, groups, 7), (2, 1, 1), (3, groups, 9), (4, 1, 1)] { + dispatch(&mut encoder, &reduction.pipelines[index], &group, x, y); + } + let owned = runtime.device().create_buffer(&wgpu::BufferDescriptor { + label: None, + size: u64::from(count) * 24 + 32, + usage: wgpu::BufferUsages::COPY_DST | wgpu::BufferUsages::COPY_SRC, + mapped_at_creation: false, + }); + encoder.copy_buffer_to_buffer(&queries, 0, &owned, 0, u64::from(count) * 16); + encoder.copy_buffer_to_buffer(&matches, 0, &owned, u64::from(count) * 16, u64::from(count) * 8); + encoder.copy_buffer_to_buffer(&selection, 0, &owned, u64::from(count) * 24, 32); + scratch.submit_and_recycle(encoder); + let words = read(runtime, &owned, owned.size()); + let query_words = &words[..count as usize * 4]; + let queries = query_words + .chunks_exact(4) + .map(|p| std::array::from_fn(|i| f32::from_bits(p[i]))) + .collect(); + let matches = words[count as usize * 4..count as usize * 6] + .chunks_exact(2) + .map(|p| [p[0], p[1]]) + .collect(); + let selection = words[count as usize * 6..].to_vec(); + let raw = read(runtime, result.buffer(), 176); + if precision == GpuIcpTrimPrecision::DoubleSingle { + assert!(raw[13..].chunks_exact(2).any(|p| p[0] != 0)); + } + let statistics = result.download(runtime).unwrap(); + assert_eq!( + statistics, + input + .trimmed_moments_resident_with_precision(runtime, pose, 0.25, 0.8, precision) + .unwrap() + .download(runtime) + .unwrap() + ); + (statistics, queries, matches, selection) +} + +#[test] +fn precise_moments_match_fixed_pose_independent_f64_oracles() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let n = 262145usize; + let target: Vec<_> = (0..n + 17) + .map(|i| [(i % 32) as f32 * 0.5, (i / 32 % 32) as f32 * 0.5, (i / 1024) as f32 * 0.5]) + .collect(); + let source: Vec<_> = target[..n] + .iter() + .enumerate() + .map(|(i, p)| { + [ + p[0] + if i % 7 == 0 { 0.1875 } else { 0.03125 }, + p[1] - 0.015625, + p[2] + 0.0078125 + if i % 11 == 0 { 0.375 } else { 0.0 }, + ] + }) + .collect(); + let target_tree = tree(&target); + let input = GpuIcpInput::upload(&runtime, &source, &target_tree).unwrap(); + let poses = [ + ("identity", Isometry3::identity()), + ( + "cpu1", + Isometry3::new( + Quat::new(-2.984037e-8, -2.9834867e-11, 1.7670175e-12, 1.0), + Vec3::new(-0.031249996, 0.015621266, -0.004607767), + ), + ), + ( + "cpu2", + Isometry3::new( + Quat::new(-9.731036e-7, -6.82904e-10, -5.5494216e-13, 1.0), + Vec3::new(-0.031249914, 0.015503033, -0.0045932755), + ), + ), + ]; + let mut maximum_mean_error = 0.0_f64; + let mut maximum_covariance_error = 0.0_f64; + for (label, pose) in poses { + let mut previous = None; + for precision in [GpuIcpTrimPrecision::CompensatedF32, GpuIcpTrimPrecision::DoubleSingle] { + let (actual, queries, matches, selection) = snapshot(&runtime, &input, pose, precision); + let mut gpu_pairs: Vec<_> = matches + .iter() + .enumerate() + .filter(|(_, m)| m[0] != u32::MAX) + .map(|(i, m)| (m[1] & 0x7fffffff, i, m[0] as usize)) + .collect(); + gpu_pairs.sort_unstable(); + let k = (gpu_pairs.len() as f64 * 0.8).floor() as usize; + assert_eq!(selection[1] as usize, k); + assert_eq!( + (selection[3], selection[4] as usize), + (gpu_pairs[k - 1].0, gpu_pairs[k - 1].1) + ); + let mean_s: [f64; 3] = std::array::from_fn(|a| { + gpu_pairs[..k].iter().map(|p| f64::from(queries[p.1][a])).sum::() / k as f64 + }); + let mean_t: [f64; 3] = std::array::from_fn(|a| { + gpu_pairs[..k].iter().map(|p| f64::from(target[p.2][a])).sum::() / k as f64 + }); + let h: [[f64; 3]; 3] = std::array::from_fn(|a| { + std::array::from_fn(|b| { + gpu_pairs[..k] + .iter() + .map(|p| { + (f64::from(queries[p.1][a]) - mean_s[a]) + * (f64::from(target[p.2][b]) - mean_t[b]) + }) + .sum() + }) + }); + let mean_error = (0..3) + .flat_map(|a| { + [ + (actual.moments.source_centroid[a] - mean_s[a]).abs(), + (actual.moments.target_centroid[a] - mean_t[a]).abs(), + ] + }) + .fold(0.0, f64::max); + let covariance_error = (0..3) + .flat_map(|a| { + (0..3).map(move |b| (actual.moments.cross_covariance[a][b] - h[a][b]).abs()) + }) + .fold(0.0, f64::max); + if let Some((old_queries, old_matches, old_selection, old_error)) = &previous { + assert_eq!(&queries, old_queries); + assert_eq!(&matches, old_matches); + assert_eq!(&selection, old_selection); + assert_eq!(&actual.moments.squared_error_sum, old_error); + } else { + previous = Some(( + queries.clone(), + matches.clone(), + selection.clone(), + actual.moments.squared_error_sum, + )); + } + if precision == GpuIcpTrimPrecision::DoubleSingle { + assert!(mean_error < 1e-9, "label={label} mean_error={mean_error}"); + for a in 0..3 { + for b in 0..3 { + let absolute = gpu_pairs[..k] + .iter() + .map(|p| { + ((f64::from(queries[p.1][a]) - mean_s[a]) + * (f64::from(target[p.2][b]) - mean_t[b])) + .abs() + }) + .sum::(); + assert!( + (actual.moments.cross_covariance[a][b] - h[a][b]).abs() + < 1e-9 + 5e-12 * absolute, + "label={label} a={a} b={b}" + ); + } + } + maximum_mean_error = maximum_mean_error.max(mean_error); + maximum_covariance_error = maximum_covariance_error.max(covariance_error); + } + } + } + eprintln!("GPU_ICP_PRECISE_ORACLES poses=3 source_points=262145 maximum_mean_error={maximum_mean_error} maximum_covariance_error={maximum_covariance_error}"); +} diff --git a/crates/spatialrust-gpu/src/kernels/icp_trim_tests.rs b/crates/spatialrust-gpu/src/kernels/icp_trim_tests.rs index 95e1c6aa..90583cd0 100644 --- a/crates/spatialrust-gpu/src/kernels/icp_trim_tests.rs +++ b/crates/spatialrust-gpu/src/kernels/icp_trim_tests.rs @@ -1,6 +1,12 @@ use super::*; use spatialrust_math::{Quat, Vec3}; +#[path = "icp_trim_precision_tests.rs"] +mod precision_tests; + +#[path = "icp_trim_math_tests.rs"] +mod math_tests; + fn read(runtime: &WgpuRuntime, buffer: &wgpu::Buffer, bytes: u64) -> Vec { let staging = runtime.device().create_buffer(&wgpu::BufferDescriptor { label: Some("trim-test-readback"), @@ -515,3 +521,11 @@ fn trimmed_gate_overflow_and_malformed_fixed_readbacks_reject_all_payload() { assert!(result.download(&runtime).is_err(), "malformed case {case}"); } } + +#[test] +fn explicit_f32_shader_preserves_the_pinned_previous_source() { + assert_eq!( + trimmed_shader(), + include_str!("../../tests/fixtures/icp_trimmed_f32_previous.wgsl").replace("\r\n", "\n") + ); +} diff --git a/crates/spatialrust-gpu/src/kernels/mod.rs b/crates/spatialrust-gpu/src/kernels/mod.rs index e7052ca9..7b85e064 100644 --- a/crates/spatialrust-gpu/src/kernels/mod.rs +++ b/crates/spatialrust-gpu/src/kernels/mod.rs @@ -60,6 +60,7 @@ pub(crate) mod icp_reduction; pub use gicp_reduction::{GpuGicpEquations, GpuGicpInput, GpuGicpNormalEquations}; pub(crate) mod iss; pub(crate) mod mls; +pub use icp_reduction::trim::GpuIcpTrimPrecision; pub use icp_reduction::{ GpuIcpInput, GpuIcpMoments, GpuIcpStatistics, GpuIcpTrimmedMoments, GpuIcpTrimmedStatistics, }; diff --git a/crates/spatialrust-gpu/src/lib.rs b/crates/spatialrust-gpu/src/lib.rs index cca83cad..4a777cfe 100644 --- a/crates/spatialrust-gpu/src/lib.rs +++ b/crates/spatialrust-gpu/src/lib.rs @@ -53,12 +53,12 @@ pub use device::{GpuDevice, WgpuDevice}; #[cfg(feature = "gpu-wgpu")] pub use kernels::{ GpuFpfhDescriptors, GpuFpfhInput, GpuFpfhNeighborhoods, GpuGicpEquations, GpuGicpInput, - GpuGicpNormalEquations, GpuIcpInput, GpuIcpMoments, GpuIcpStatistics, GpuIcpTrimmedMoments, - GpuIcpTrimmedStatistics, GpuIssInput, GpuIssKeypoints, GpuIssNeighborhoods, GpuMlsDiagnostics, - GpuMlsInput, GpuMlsMode, GpuMlsNeighborhoods, GpuMlsPositions, GpuMlsReadback, - GpuNeighborhoodKind, GpuNeighborhoodReadback, GpuNeighborhoods, GpuPointToPlaneEquations, - GpuPointToPlaneInput, GpuPointToPlaneNormalEquations, GpuQueryBatch, GPU_FPFH_DESCRIPTOR_LEN, - MAX_GPU_K_NEIGHBORS, + GpuGicpNormalEquations, GpuIcpInput, GpuIcpMoments, GpuIcpStatistics, GpuIcpTrimPrecision, + GpuIcpTrimmedMoments, GpuIcpTrimmedStatistics, GpuIssInput, GpuIssKeypoints, + GpuIssNeighborhoods, GpuMlsDiagnostics, GpuMlsInput, GpuMlsMode, GpuMlsNeighborhoods, + GpuMlsPositions, GpuMlsReadback, GpuNeighborhoodKind, GpuNeighborhoodReadback, + GpuNeighborhoods, GpuPointToPlaneEquations, GpuPointToPlaneInput, + GpuPointToPlaneNormalEquations, GpuQueryBatch, GPU_FPFH_DESCRIPTOR_LEN, MAX_GPU_K_NEIGHBORS, }; #[cfg(feature = "gpu-wgpu")] diff --git a/crates/spatialrust-gpu/src/shaders/icp_trim_moments.wgsl b/crates/spatialrust-gpu/src/shaders/icp_trim_moments.wgsl new file mode 100644 index 00000000..ccf1d0e1 --- /dev/null +++ b/crates/spatialrust-gpu/src/shaders/icp_trim_moments.wgsl @@ -0,0 +1,63 @@ +// Carry two f32 words through normalization, mean division and products. +// Matching/queries stay f32. Existing add/merge retain compensated reductions. +// Host always supplies params.unused=0. The runtime integer operation preserves +// each rounded f32 value while preventing elimination of error-free transforms +// through floating-point reassociation across their dependent expressions. +fn ds_round(value: f32) -> f32 { + return bitcast(bitcast(value)^params.unused); +} +fn precise_add(a: vec2, b: f32) -> vec2 { + let s=ds_round(a.x+b); let v=ds_round(s-a.x); + let e=ds_round(ds_round(a.x-ds_round(s-v))+ds_round(b-v))+a.y; + let h=ds_round(s+e); + return vec2(h,ds_round(e-ds_round(h-s))); +} +fn precise_merge(a: vec2, b: vec2) -> vec2 { + return precise_add(precise_add(a,b.x),b.y); +} +fn ds_multiply_words(a: u32, b: u32) -> vec2 { + let a0=a&65535u; let a1=a>>16u; let b0=b&65535u; let b1=b>>16u; + let w0=a0*b0; let t=a1*b0+(w0>>16u); let w1=a0*b1+(t&65535u); + return vec2((w1<<16u)|(w0&65535u),a1*b1+(t>>16u)+(w1>>16u)); +} +fn ds_scale_word(mantissa: u32, exponent: i32) -> f32 { + if mantissa==0u { return 0.0; } + let bits=firstLeadingBit(mantissa)+1u; + let fraction=f32(mantissa)*bitcast((127u-bits)<<23u); + let combined=exponent+i32(bits); + if combined>128 { return bitcast(0x7f800000u); } + return ldexp(fraction,combined); +} +// Integer 24x24 multiplication prevents floating-point contraction or an +// unfused WGSL fma from erasing the exact product's low word. +fn two_product(a: f32, b: f32) -> vec2 { + let x=bitcast(a); let y=bitcast(b); + let xe=(x>>23u)&255u; let ye=(y>>23u)&255u; + if xe==255u || ye==255u { return vec2(bitcast(0x7f800000u),0.0); } + let xm=(x&0x7fffffu)|select(0x800000u,0u,xe==0u); + let ym=(y&0x7fffffu)|select(0x800000u,0u,ye==0u); + let exponent=select(i32(xe)-150,-149,xe==0u)+select(i32(ye)-150,-149,ye==0u); + let product=ds_multiply_words(xm,ym); + let upper=(product.x>>24u)|(product.y<<8u); + let lower=product.x&0xffffffu; + let sign=(x^y)&0x80000000u; + let hi=bitcast(bitcast(ds_scale_word(upper,exponent+24))|sign); + let lo=bitcast(bitcast(ds_scale_word(lower,exponent))|sign); + return vec2(hi,lo); +} +fn ds_divide(a: vec2, divisor: f32) -> vec2 { + let quotient=a.x/divisor; + let product=two_product(quotient,divisor); + let remainder=precise_merge(a,-product); + let correction=(remainder.x+remainder.y)/divisor; + return precise_add(vec2(quotient,0.0),correction); +} +fn ds_offset(point: f32, anchor: f32, scale: f32) -> vec2 { + return ds_divide(precise_add(vec2(point,0.0),-anchor),scale); +} +fn ds_product(a: vec2, b: vec2) -> vec2 { + var total=two_product(a.x,b.x); + total=precise_add(total,a.x*b.y); + total=precise_add(total,a.y*b.x); + return precise_add(total,a.y*b.y); +} diff --git a/crates/spatialrust-gpu/tests/fixtures/icp_trimmed_f32_previous.wgsl b/crates/spatialrust-gpu/tests/fixtures/icp_trimmed_f32_previous.wgsl new file mode 100644 index 00000000..ae5b07ae --- /dev/null +++ b/crates/spatialrust-gpu/tests/fixtures/icp_trimmed_f32_previous.wgsl @@ -0,0 +1,145 @@ +struct Params { + row0: vec4, row1: vec4, row2: vec4, + source_anchor: vec4, target_anchor: vec4, + count: u32, targets: u32, groups: u32, unused: u32, + bound: f32, max_distance: f32, scale: f32, padding: f32, +} +struct Neighbor { index: u32, distance: f32 } +struct Partial { values: array,16>, count: u32, pad0: u32, pad1: u32, pad2: u32 } +struct Statistics { + status: atomic, count: u32, reserved0: u32, reserved1: u32, + source_anchor: vec4, target_anchor: vec4, values: array,16>, +} +@group(0) @binding(0) var params: Params; +@group(0) @binding(1) var source: array>; +@group(0) @binding(2) var target_points: array>; +@group(0) @binding(3) var queries: array>; +@group(0) @binding(4) var matches: array; +@group(0) @binding(5) var partials: array; +@group(0) @binding(6) var result: Statistics; +var sums: array,64>; +var counts: array; + +fn finite(v: f32) -> bool { return abs(v)<=3.4028234663852886e38; } +fn finite3(v: vec3) -> bool { return finite(v.x) && finite(v.y) && finite(v.z); } +fn dot3(a: vec3, b: vec3) -> f32 { return a.x*b.x+a.y*b.y+a.z*b.z; } +fn transform_point(p: vec3) -> vec3 { + return vec3(dot3(params.row0.xyz,p)+params.row0.w, + dot3(params.row1.xyz,p)+params.row1.w,dot3(params.row2.xyz,p)+params.row2.w); +} +fn add(a: vec2, b: f32) -> vec2 { + let s=a.x+b; let v=s-a.x; let e=(a.x-(s-v))+(b-v)+a.y; let h=s+e; + return vec2(h,e-(h-s)); +} +fn merge(a: vec2, b: vec2) -> vec2 { return add(add(a,b.x),b.y); } +fn valid_match(i: u32) -> bool { + let m=matches[i]; + if m.index==0xffffffffu { + if m.distance!=params.bound { atomicOr(&result.status,2u); } + return false; + } + if m.index>=params.targets || !(m.distance>=0.0 && m.distance<=params.bound) { + atomicOr(&result.status,2u); return false; + } + return true; +} + +@compute @workgroup_size(64) +fn transform_source(@builtin(global_invocation_id) id: vec3) { + if id.x>=params.count { return; } + let p=transform_point(source[id.x].xyz); + if !finite3(p) { atomicOr(&result.status,1u); queries[id.x]=vec4(0.0); } + else { queries[id.x]=vec4(p,0.0); } +} + +@compute @workgroup_size(64) +fn centroid_partials(@builtin(workgroup_id) group: vec3, @builtin(local_invocation_index) lane: u32) { + let term=select(group.y,15u,group.y==6u); + let anchor=transform_point(params.source_anchor.xyz); + var total=vec2(0.0); var count=0u; + for (var i=group.x*64u+lane;i0u;stride/=2u) { + if lane) { + if id.x!=0u { return; } + var count=0u; + for (var g=0u;g(transform_point(params.source_anchor.xyz),0.0); + result.target_anchor=params.target_anchor; + if !finite3(result.source_anchor.xyz) { atomicOr(&result.status,1u); } + for (var term=0u;term<16u;term+=1u) { + result.values[term]=vec2(0.0); + if term<6u || term==15u { + var total=vec2(0.0); + for (var g=0u;g0u || term==15u { + if term<6u { result.values[term]=merge(vec2(total.x/f32(count),0.0),vec2(total.y/f32(count),0.0)); } + else { result.values[term]=total; } + } + } + } +} + +@compute @workgroup_size(64) +fn covariance_partials(@builtin(workgroup_id) group: vec3, @builtin(local_invocation_index) lane: u32) { + let term=group.y; let row=term/3u; let col=term%3u; + let mean_s=result.values[row].x+result.values[row].y; + let mean_t=result.values[col+3u].x+result.values[col+3u].y; + var total=vec2(0.0); + if result.count>0u && atomicLoad(&result.status)==0u { + for (var i=group.x*64u+lane;i0u;stride/=2u) { + if lane) { + if id.x!=0u { return; } + for (var term=6u;term<15u;term+=1u) { + var total=vec2(0.0); + for (var g=0u;g(0.0); result.target_anchor=vec4(0.0); } + if atomicLoad(&result.status)!=0u || result.reserved0==0u { + result.count=0u; result.reserved0=0u; result.source_anchor=vec4(0.0); result.target_anchor=vec4(0.0); + for (var term=0u;term<16u;term+=1u) { result.values[term]=vec2(0.0); } + } +} + +struct Selection { gated:u32, retained:u32, rank:u32, distance:u32, source_index:u32, unused0:u32, unused1:u32, unused2:u32 } +@group(0) @binding(7) var selection: Selection; +fn retained_match(i:u32)->bool { let distance=bitcast(matches[i].distance)&0x7fffffffu; return selection.retained>0u && (distance SpatialResult { - self.align_reduction_impl(runtime, source, target, false) + self.align_reduction_impl(runtime, source, target, false, GpuIcpTrimPrecision::DoubleSingle) } /// Aligns using GPU stable distance/source-index trimming and resident @@ -54,7 +54,28 @@ impl IcpRegistration { source: &PointCloud, target: &PointCloud, ) -> SpatialResult { - self.align_reduction_impl(runtime, source, target, true) + self.align_with_gpu_trimmed_reduction_with_precision( + runtime, + source, + target, + GpuIcpTrimPrecision::DoubleSingle, + ) + } + + /// Aligns with an explicit retained-moment precision choice. DoubleSingle + /// (the default resident trimmed route) carries normalized quotient/product + /// residuals in two f32 words. CompensatedF32 retains the earlier lower-cost + /// arithmetic. Transform/search precision, all-gated fitness, CPU solving + /// and explicit fixed-size transfer receipts are the same for both choices. + #[cfg(feature = "register-icp-gpu")] + pub fn align_with_gpu_trimmed_reduction_with_precision( + &self, + runtime: &WgpuRuntime, + source: &PointCloud, + target: &PointCloud, + precision: GpuIcpTrimPrecision, + ) -> SpatialResult { + self.align_reduction_impl(runtime, source, target, true, precision) } fn align_reduction_impl( @@ -63,6 +84,7 @@ impl IcpRegistration { source: &PointCloud, target: &PointCloud, trimmed: bool, + precision: GpuIcpTrimPrecision, ) -> SpatialResult { self.validate_inputs(source, target)?; if !trimmed && self.trim_fraction != 1.0 { @@ -77,11 +99,12 @@ impl IcpRegistration { let mut evaluation_batches = 0; let mut evaluate = |pose| -> SpatialResult { let statistics = input - .trimmed_moments_resident( + .trimmed_moments_resident_with_precision( runtime, pose, self.config.max_correspondence_distance, self.trim_fraction, + precision, )? .download(runtime)?; transfers.record( diff --git a/crates/spatialrust-registration/src/lib.rs b/crates/spatialrust-registration/src/lib.rs index 1237da83..a64ce734 100644 --- a/crates/spatialrust-registration/src/lib.rs +++ b/crates/spatialrust-registration/src/lib.rs @@ -51,6 +51,9 @@ pub use fpfh::{ #[cfg(feature = "register-icp-gpu")] pub use icp::{GpuIcpDiagnostics, GpuIcpReductionDiagnostics}; +#[cfg(feature = "register-icp-gpu")] +pub use spatialrust_gpu::GpuIcpTrimPrecision; + #[cfg(feature = "register-gicp-gpu")] pub use gicp::{GpuGicpDiagnostics, GpuGicpIteration, GpuGicpReductionDiagnostics}; diff --git a/crates/spatialrust-registration/tests/gpu_icp_trimmed_precision.rs b/crates/spatialrust-registration/tests/gpu_icp_trimmed_precision.rs new file mode 100644 index 00000000..655ce061 --- /dev/null +++ b/crates/spatialrust-registration/tests/gpu_icp_trimmed_precision.rs @@ -0,0 +1,153 @@ +#![cfg(all(feature = "register-icp-gpu", not(target_arch = "wasm32")))] +#[path = "gpu_icp_trimmed_reduction.rs"] +mod reference; +use spatialrust_gpu::WgpuRuntime; +use spatialrust_registration::{ + GpuIcpTrimPrecision, IcpConfig, IcpConvergenceCriteria, IcpRegistration, +}; + +pub fn registration(iterations: usize, fraction: f64) -> IcpRegistration { + IcpRegistration::new(IcpConfig { + max_iterations: iterations, + max_correspondence_distance: 0.25, + ..Default::default() + }) + .with_trim_fraction(fraction) + .with_convergence_criteria(IcpConvergenceCriteria { + translation_epsilon: 0.0, + rotation_epsilon: 0.0, + fitness_epsilon: 0.0, + }) +} + +#[test] +#[ignore = "large native GPU precision regressions; run explicitly in release mode"] +fn precision_large_clouds_preserve_pose_counts_fitness_and_iteration_trajectories() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let (mut maximum, mut previous_maximum, mut conditions) = (0.0_f32, 0.0_f32, 0); + for n in [4097, 65537, 262145] { + for profile in 0..3 { + let fraction = [0.5, 0.8, 0.9][profile]; + let source = reference::fixture(n, profile, 1.0, true); + let target = reference::fixture(n + 17, profile, 1.0, false); + for iterations in 1..=3 { + let registration = registration(iterations, fraction); + let cpu = registration.align_with_trace(&source, &target).unwrap(); + let precise = registration + .align_with_gpu_trimmed_reduction(&runtime, &source, &target) + .unwrap(); + let previous = registration + .align_with_gpu_trimmed_reduction_with_precision( + &runtime, + &source, + &target, + GpuIcpTrimPrecision::CompensatedF32, + ) + .unwrap(); + let error = reference::pose_error( + &source, + cpu.result.transform, + precise.diagnostics.result.transform, + ); + maximum = maximum.max(error); + previous_maximum = previous_maximum.max(reference::pose_error( + &source, + cpu.result.transform, + previous.diagnostics.result.transform, + )); + assert!( + error < 3e-5, + "n={n} profile={profile} iterations={iterations} error={error}" + ); + assert_eq!(precise.diagnostics.stop_reason, cpu.stop_reason); + assert_eq!(precise.diagnostics.result.converged, cpu.result.converged); + for (a, b) in precise.diagnostics.history.iter().zip(&cpu.history) { + assert_eq!( + (a.correspondences, a.evaluated_correspondences), + (b.correspondences, b.evaluated_correspondences) + ); + } + let (count, fitness) = reference::independent_fitness( + &source, + &target, + precise.diagnostics.result.transform, + 0.25, + ); + assert_eq!( + count, + precise.diagnostics.history.last().unwrap().evaluated_correspondences + ); + assert!( + (fitness - precise.diagnostics.result.fitness).abs() <= 1e-10 + fitness * 1e-4 + ); + assert_eq!(precise.evaluation_batches, iterations + 1); + assert_eq!(precise.transfers.device_to_host_bytes(), 176 * (iterations + 1) as u64); + assert_eq!(precise.transfers, previous.transfers); + conditions += 1; + } + } + } + eprintln!("GPU_ICP_PRECISE_TRAJECTORIES conditions={conditions} maximum_pose_error={maximum} previous_maximum_pose_error={previous_maximum}"); +} + +#[test] +#[ignore = "requires a native GPU adapter; run explicitly"] +fn precision_choices_keep_fraction_one_and_mixed_concurrent_results() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let source = reference::fixture(4097, 1, 1.0, true); + let target = reference::fixture(4114, 1, 1.0, false); + let untrimmed = registration(3, 1.0); + let ordinary = untrimmed.align_with_gpu_reduction(&runtime, &source, &target).unwrap(); + for precision in [GpuIcpTrimPrecision::CompensatedF32, GpuIcpTrimPrecision::DoubleSingle] { + assert_eq!( + untrimmed + .align_with_gpu_trimmed_reduction_with_precision( + &runtime, &source, &target, precision + ) + .unwrap(), + ordinary + ); + } + let registration = registration(3, 0.8); + let expected = + [GpuIcpTrimPrecision::CompensatedF32, GpuIcpTrimPrecision::DoubleSingle].map(|p| { + registration + .align_with_gpu_trimmed_reduction_with_precision(&runtime, &source, &target, p) + .unwrap() + }); + assert_eq!( + expected[1], + registration.align_with_gpu_trimmed_reduction(&runtime, &source, &target).unwrap() + ); + std::thread::scope(|scope| { + let workers: Vec<_> = (0..4) + .map(|i| { + let registration = ®istration; + let runtime = &runtime; + let source = &source; + let target = ⌖ + let expected = &expected; + scope.spawn(move || { + for round in 0..2 { + let index = (i + round) % 2; + let precision = [ + GpuIcpTrimPrecision::CompensatedF32, + GpuIcpTrimPrecision::DoubleSingle, + ][index]; + assert_eq!( + registration + .align_with_gpu_trimmed_reduction_with_precision( + runtime, source, target, precision + ) + .unwrap(), + expected[index] + ); + } + }) + }) + .collect(); + for worker in workers { + worker.join().unwrap(); + } + }); +} diff --git a/crates/spatialrust-registration/tests/gpu_icp_trimmed_precision_benchmark.rs b/crates/spatialrust-registration/tests/gpu_icp_trimmed_precision_benchmark.rs new file mode 100644 index 00000000..da93ac26 --- /dev/null +++ b/crates/spatialrust-registration/tests/gpu_icp_trimmed_precision_benchmark.rs @@ -0,0 +1,202 @@ +#![cfg(all(feature = "register-icp-gpu", not(target_arch = "wasm32")))] +#[path = "gpu_icp_trimmed_reduction.rs"] +mod reference; +use spatialrust_core::HasPositions3; +use spatialrust_gpu::{GpuIcpInput, WgpuRuntime}; +use spatialrust_registration::GpuIcpTrimPrecision; +use spatialrust_search::KdTree; +use std::time::Instant; + +#[test] +#[ignore = "isolated completed-call precision survey; run this exact test alone in release mode"] +fn completed_trimmed_icp_precision_calls() { + let runtime = WgpuRuntime::new_headless().unwrap(); + let info = runtime.adapter_info(); + println!("GPU_ICP_PRECISE_ENV backend={} device_type={}", info.backend, info.device_type); + for n in [257, 4097, 65537, 262145] { + for profile in 0..3 { + let fraction = [0.5, 0.8, 0.9][profile]; + let source = reference::fixture(n, profile, 1.0, true); + let target = reference::fixture(n + 17, profile, 1.0, false); + let registration = reference::registration(1.0, fraction, 0); + let cpu = registration.align_with_trace(&source, &target).unwrap(); + let legacy = + registration.align_with_gpu_correspondences(&runtime, &source, &target).unwrap(); + let fast = registration + .align_with_gpu_trimmed_reduction_with_precision( + &runtime, + &source, + &target, + GpuIcpTrimPrecision::CompensatedF32, + ) + .unwrap(); + let precise = + registration.align_with_gpu_trimmed_reduction(&runtime, &source, &target).unwrap(); + let fast_error = reference::pose_error( + &source, + cpu.result.transform, + fast.diagnostics.result.transform, + ); + let precise_error = reference::pose_error( + &source, + cpu.result.transform, + precise.diagnostics.result.transform, + ); + assert!(precise_error < 3e-5, "n={n} profile={profile} error={precise_error}"); + for diagnostics in [&legacy.diagnostics, &fast.diagnostics, &precise.diagnostics] { + assert_eq!(diagnostics.stop_reason, cpu.stop_reason); + assert_eq!(diagnostics.result.iterations, 3); + assert_eq!(diagnostics.result.converged, cpu.result.converged); + for (a, b) in diagnostics.history.iter().zip(&cpu.history) { + assert_eq!( + (a.correspondences, a.evaluated_correspondences), + (b.correspondences, b.evaluated_correspondences) + ); + } + let (count, fitness) = reference::independent_fitness( + &source, + &target, + diagnostics.result.transform, + 0.25, + ); + assert_eq!(count, diagnostics.history.last().unwrap().evaluated_correspondences); + assert!((fitness - diagnostics.result.fitness).abs() <= 1e-10 + fitness * 1e-4); + } + assert_eq!(fast.transfers, precise.transfers); + assert_eq!(precise.evaluation_batches, 4); + let (sx, sy, sz) = source.positions3().unwrap(); + let (tx, ty, tz) = target.positions3().unwrap(); + let positions: Vec<_> = (0..n).map(|i| [sx[i], sy[i], sz[i]]).collect(); + let tree = KdTree::from_slices(tx, ty, tz); + let input = GpuIcpInput::upload(&runtime, &positions, &tree).unwrap(); + let pose = precise.diagnostics.result.transform; + let prepared = [GpuIcpTrimPrecision::CompensatedF32, GpuIcpTrimPrecision::DoubleSingle] + .map(|p| { + input + .trimmed_moments_resident_with_precision(&runtime, pose, 0.25, fraction, p) + .unwrap() + .download(&runtime) + .unwrap() + }); + assert_eq!( + prepared[0].evaluated_correspondences, + prepared[1].evaluated_correspondences + ); + assert_eq!(prepared[0].moments.correspondences, prepared[1].moments.correspondences); + assert_eq!( + prepared[0].moments.squared_error_sum, + prepared[1].moments.squared_error_sum + ); + for _ in 0..3 { + assert_eq!(registration.align_with_trace(&source, &target).unwrap(), cpu); + assert_eq!( + registration + .align_with_gpu_correspondences(&runtime, &source, &target) + .unwrap(), + legacy + ); + assert_eq!( + registration + .align_with_gpu_trimmed_reduction_with_precision( + &runtime, + &source, + &target, + GpuIcpTrimPrecision::CompensatedF32 + ) + .unwrap(), + fast + ); + assert_eq!( + registration + .align_with_gpu_trimmed_reduction(&runtime, &source, &target) + .unwrap(), + precise + ); + for (index, precision) in + [GpuIcpTrimPrecision::CompensatedF32, GpuIcpTrimPrecision::DoubleSingle] + .into_iter() + .enumerate() + { + assert_eq!( + input + .trimmed_moments_resident_with_precision( + &runtime, pose, 0.25, fraction, precision + ) + .unwrap() + .download(&runtime) + .unwrap(), + prepared[index] + ); + } + } + let ( + mut cpu_ms, + mut legacy_ms, + mut fast_ms, + mut precise_ms, + mut prepared_fast_ms, + mut prepared_precise_ms, + ) = (vec![], vec![], vec![], vec![], vec![], vec![]); + for round in 0..20 { + for offset in 0..6 { + let stage = (round + offset) % 6; + let start = Instant::now(); + match stage { + 0 => { + let output = registration.align_with_trace(&source, &target).unwrap(); + cpu_ms.push(start.elapsed().as_secs_f64() * 1000.0); + assert_eq!(output, cpu); + } + 1 => { + let output = registration + .align_with_gpu_correspondences(&runtime, &source, &target) + .unwrap(); + legacy_ms.push(start.elapsed().as_secs_f64() * 1000.0); + assert_eq!(output, legacy); + } + 2 => { + let output = registration + .align_with_gpu_trimmed_reduction_with_precision( + &runtime, + &source, + &target, + GpuIcpTrimPrecision::CompensatedF32, + ) + .unwrap(); + fast_ms.push(start.elapsed().as_secs_f64() * 1000.0); + assert_eq!(output, fast); + } + 3 => { + let output = registration + .align_with_gpu_trimmed_reduction(&runtime, &source, &target) + .unwrap(); + precise_ms.push(start.elapsed().as_secs_f64() * 1000.0); + assert_eq!(output, precise); + } + _ => { + let index = stage - 4; + let precision = [ + GpuIcpTrimPrecision::CompensatedF32, + GpuIcpTrimPrecision::DoubleSingle, + ][index]; + let output = input + .trimmed_moments_resident_with_precision( + &runtime, pose, 0.25, fraction, precision, + ) + .unwrap(); + runtime.wait_idle(); + let elapsed = start.elapsed().as_secs_f64() * 1000.0; + if index == 0 { + prepared_fast_ms.push(elapsed); + } else { + prepared_precise_ms.push(elapsed); + } + assert_eq!(output.download(&runtime).unwrap(), prepared[index]); + } + } + } + } + println!("GPU_ICP_PRECISE_BENCH n={n} profile={profile} fraction={fraction} iterations=3 fast_pose_error={fast_error} precise_pose_error={precise_error} cpu_fitness={} fast_fitness={} precise_fitness={} retained={} gated={} legacy_h2d={} legacy_d2h={} precise_h2d={} precise_d2h={} cpu_ms={cpu_ms:?} legacy_ms={legacy_ms:?} fast_ms={fast_ms:?} precise_ms={precise_ms:?} prepared_fast_ms={prepared_fast_ms:?} prepared_precise_ms={prepared_precise_ms:?}",cpu.result.fitness,fast.diagnostics.result.fitness,precise.diagnostics.result.fitness,prepared[1].moments.correspondences,prepared[1].evaluated_correspondences,legacy.transfers.host_to_device_bytes(),legacy.transfers.device_to_host_bytes(),precise.transfers.host_to_device_bytes(),precise.transfers.device_to_host_bytes()); + } + } +} diff --git a/crates/spatialrust-registration/tests/gpu_icp_trimmed_reduction_benchmark.rs b/crates/spatialrust-registration/tests/gpu_icp_trimmed_reduction_benchmark.rs index f5cfe959..6cc27469 100644 --- a/crates/spatialrust-registration/tests/gpu_icp_trimmed_reduction_benchmark.rs +++ b/crates/spatialrust-registration/tests/gpu_icp_trimmed_reduction_benchmark.rs @@ -60,11 +60,7 @@ fn completed_resident_trimmed_icp_calls() { ); assert!(maximum.is_finite()); let within_pose_tolerance = maximum < 3e-5; - // Preserve adverse large gated cases in the measurements rather - // than silently increasing the small-fixture agreement tolerance. - if n != 262145 || profile != 1 { - assert!(within_pose_tolerance, "n={n} profile={profile} max_pose_error={maximum}"); - } + assert!(within_pose_tolerance, "n={n} profile={profile} max_pose_error={maximum}"); assert_eq!(expected_gpu.diagnostics.result.iterations, 3); assert_eq!(expected_gpu.evaluation_batches, 4); for (a, b) in expected_cpu.history.iter().zip(&expected_gpu.diagnostics.history) { diff --git a/crates/spatialrust/src/lib.rs b/crates/spatialrust/src/lib.rs index 561992ac..ac7f47f2 100644 --- a/crates/spatialrust/src/lib.rs +++ b/crates/spatialrust/src/lib.rs @@ -242,7 +242,7 @@ pub use spatialrust_registration::{PointToPlaneIcp, PointToPlaneIcpConfig}; pub use spatialrust_registration::{GpuPointToPlaneIteration, GpuPointToPlaneReductionDiagnostics}; #[cfg(feature = "register-icp-gpu")] -pub use spatialrust_registration::GpuIcpReductionDiagnostics; +pub use spatialrust_registration::{GpuIcpReductionDiagnostics, GpuIcpTrimPrecision}; #[cfg(feature = "register-gicp")] pub use spatialrust_registration::{GicpConfig, GicpRegistration}; diff --git a/docs/ARCHITECTURE.md b/docs/ARCHITECTURE.md index 769e642d..35969e17 100644 --- a/docs/ARCHITECTURE.md +++ b/docs/ARCHITECTURE.md @@ -42,10 +42,22 @@ stopping denominator. A 176-byte result adds the all-gated count in reserved wor Private histogram/selection scratch stays leased through the same submission. CPU target-tree construction/Horn solving/composition/convergence remain; there is no payload readback or fallback. GPU distance rounding can change selection -near ties compared with CPU, even when retained counts agree. A measured large -partially gated case exceeds CPU pose-agreement tolerance; the evidence retains -its error/fitness and does not classify its speed as accuracy-qualified. See -[resident trimming evidence](../notes/2026-10-13_gpu_icp_resident_trimming.md). +near ties compared with CPU, even when retained counts agree. See +[initial resident trimming evidence](../notes/2026-10-13_gpu_icp_resident_trimming.md). + +Resident trimming now defaults to `GpuIcpTrimPrecision::DoubleSingle`. Integer +24x24-bit products retain residuals, and a private runtime zero word in dependent +bitcasts preserves rounded f32 values through compensated arithmetic. Two f32 +words carry coordinate differences, normalization, mean division and centered +products, without relying on a fused WGSL `fma`. Explicit `*_with_precision` +methods offer the previous `CompensatedF32` arithmetic. Reduction pipelines are +lazy per runtime and precision; the unchanged selector and scratch are shared. +No new pass, buffer layout, transfer or limit is added. Fraction one delegates +to canonical ordinary moments for either choice. Transform/search and all-gated +fitness keep their f32 arithmetic, and CPU solving/composition/convergence remain. +The earlier large partially gated case now meets the original pose tolerance; +near-gate/tied/degenerate f32 matching can still change other trajectories. See +[trimming precision evidence](../notes/2026-10-13_gpu_icp_trimmed_precision.md). GICP reuses this search through its own explicit GPU entry point, retaining unbounded nearest matching before the existing host gate and ordered CPU diff --git a/docs/ROADMAP.md b/docs/ROADMAP.md index dccad123..65863477 100644 --- a/docs/ROADMAP.md +++ b/docs/ROADMAP.md @@ -1140,4 +1140,5 @@ Dependencies and explicit transfers come before broader algorithm coverage. | GICP resident radius covariances | Complete | Complete inclusive GPU radius neighborhoods with explicit capacity, compensated central sums/eigenvectors retained for equations; sparse epsilon*I and overflow rejection | Independent brute-force/f64/spectral references, 162 pose/scale/radius/stopping cases, sparse extent beyond grid capacity, receipts/ownership/status checks and completed API comparisons; `notes/2026-10-13_gpu_gicp_resident_radius_covariances.md` | | GICP resident grid covariances | Complete | Reused original positions and canonical GPU grid gather/covariance sweeps; no source KD-tree/neighborhood rows, private grid scratch and complete-support rejection | Exact legacy covariance-word parity, 162 registration conditions, boundary/error/ownership/concurrency and transfer checks with completed-call comparisons; `notes/2026-10-13_gpu_gicp_resident_grid_covariances.md` | | ICP resident trimming | Complete | Stable GPU distance/source-index selection and exact binary64 retained-count floor; selected moments plus all-gated fitness in one submission, fixed 176-byte readback | Independent comparison-sort/f64 moment references, rounding/boundary/ownership/status checks, 486 registration conditions, completed comparisons and an explicit adverse large-gate precision case; `notes/2026-10-13_gpu_icp_resident_trimming.md` | -| Remaining registration updates | Planned | Small pose-update boundaries and adverse trimming precision | Pose/convergence/fitness references, numerical conditioning, source/target residency and end-to-end registration timing | +| ICP trimmed moment precision | Complete | Default two-word coordinate/quotient/product arithmetic with integer product residuals and preserved rounding; explicit previous f32 choice, same GPU matching/readback/CPU solve | 2,076 independent arithmetic comparisons, fixed-pose f64 moments, 27 large trajectories, mixed precision/concurrency and completed timing; `notes/2026-10-13_gpu_icp_trimmed_precision.md` | +| Remaining registration updates | Planned | Small pose-update boundaries and broader numerical conditioning | Pose/convergence/fitness references, numerical conditioning, source/target residency and end-to-end registration timing | diff --git a/notes/2026-10-13_gpu_icp_trimmed_precision.md b/notes/2026-10-13_gpu_icp_trimmed_precision.md new file mode 100644 index 00000000..2285b71c --- /dev/null +++ b/notes/2026-10-13_gpu_icp_trimmed_precision.md @@ -0,0 +1,213 @@ +# Resident trimmed ICP moment precision + +Date: 2026-10-13. Named slice: accurate two-word GPU moments for trimmed ICP, +with explicit precision selection and the existing CPU pose-solve boundary. + +## API and execution + +Resident trimmed ICP now defaults to `GpuIcpTrimPrecision::DoubleSingle`. +The explicit alternatives are +`GpuIcpInput::trimmed_moments_resident_with_precision` and +`IcpRegistration::align_with_gpu_trimmed_reduction_with_precision`; choosing +`CompensatedF32` retains the initial route's coordinate/product arithmetic. +The precision enum is exported by GPU, registration and the facade under their +existing GPU feature flags. Both methods validate and propagate errors without +CPU fallback or automatic backend placement. Ordinary CPU ICP stays unchanged. + +The new route carries two f32 words through point-minus-anchor offsets, +normalization, centroid division, centering and covariance products. Integer +24x24-bit mantissa multiplication preserves the high/low product, rather than +assuming fused multiply-add behavior. Dependent bitcasts use the private +runtime uniform zero word to preserve rounded f32 intermediate values through +compensated addition/merging. The host always initializes that word to zero. +This prevents algebraic elimination of the correction terms on the measured +backend, without exposing another user parameter or copying data to CPU. + +WGSL permits floating-point reassociation, and its `fma` can be implemented as +ordinary multiplication followed by addition. An error-free product cannot +assume IEEE fused behavior. See the primary +[WGSL floating-point rules](https://www.w3.org/TR/WGSL/#floating-point-accuracy) +and [fma specification](https://www.w3.org/TR/WGSL/#fma-builtin). + +The earlier experiment's low correction words were all zero on this backend. +The preserved-rounding route retains them. At the initial large gated pose, +source/target membership, device queries and selected indices agreed exactly +with scalar matching. The initial maximum mean error fell from +0.00000043506811664428824 to 0.00000000000005684341886080802, and maximum +covariance error from 74.13663649559021 to 0.000025928020477294922 against +independent f64 sums over identical GPU-selected points. These observations +identify moment rounding as an initial source of trajectory divergence; later +f32 transforms/search can still change computed distances or selected sets. + +Five reduction pipelines are initialized lazily per precision on one runtime; +the stable distance/index selector is shared. Source/target/tree upload, +six canonical phases plus selection digits, scratch ownership, workgroup sizes, +binding/shared-memory limits, status bits and the 176-byte output layout remain. +Fraction one uses canonical untrimmed kernels for either precision and preserves +their complete diagnostics. The all-gated fitness arithmetic remains canonical; +the new arithmetic applies to retained centroids/covariance, not its denominator. +CPU retains target-tree construction, Horn solving, composition and stopping. + +There is no additional pass, point/match payload readback, host wait, repeated +upload, dependency or unsafe boundary. For S source points, T target points and +M target nodes, H2D stays `16*S+32*T+32*M`; I updates require I+1 evaluations and +`176*(I+1)` D2H bytes. Uniform controls remain excluded from payload receipts. +Both precision caches and shared scratch support concurrent callers. + +DoubleSingle improves retained-moment accuracy, but its operations and storage +still use f32 words. Subnormal/underflow behavior, extreme dynamic range, +matching near gates/ties and degenerate geometry retain numerical limits. It +does not supply native f64 device execution or guarantee CPU-identical poses +for arbitrary inputs. No change is made to source transform/search arithmetic. + +## Validation + +Ten native GPU trimming primitive tests passed, including: + +- 2,076 independent binary64 comparisons for compensated addition, exact + two-product, corrected division and two-word products: 519 input rows with + cancellation/product residuals, signs, integer boundaries and seeded varied + exponents. Maximum relative error was 0.000000000000012316562842982828; + the enforced tolerance is 1e-11. +- Three fixed-pose f64 oracles with 262,145 source points. Owned snapshots copy + queries/matches/selection before scratch recycling, compare thresholds to + independent sorting, and accumulate selected moments in f64. Queries, + matches, selected sets and all-gated error are identical between precision + choices. Across poses, maximum precise mean error was + 0.0000000000004334310688136611 and covariance error + 0.00004225969314575195. The mean limit is 1e-9; covariance uses an absolute + floor plus 5e-12 times the sum of absolute reference products. +- A frozen complete shader from source commit + `741efaf4614f0b53651ce45c119f5fd7c5d1e50e` verifies that the explicit f32 + choice retains its earlier shader expressions. Existing fraction/count + rounding, membership/ordering, all-gated fitness, ownership/concurrency, + numeric error/overflow and malformed-readback checks also pass. + +The large native registration suite checks 27 trajectories: sizes +4097/65537/262145, unique/partially gated/duplicate profiles with fractions +0.5/0.8/0.9, and one/two/three updates. Every precise result meets the unchanged +3e-5 CPU pose-coordinate agreement threshold. Maximum difference is +0.0000076293945, compared with the explicit previous choice's 0.0008154977. +Counts, independent all-gated fitness, iteration limit/nonconvergence and full +receipts agree. Fraction one is exactly canonical, and four callers alternate +both precision choices while retaining each mode's complete expected result. + +The existing 486 size/profile/scale/fraction/stopping conditions and three +rotated initial guesses passed after the precision change. Maximum +scale-normalized CPU pose difference remains 0.0000019073486. The old benchmark +now enforces its original pose threshold on the formerly adverse large case +without an exception. Facade smoke tests, all-feature registration library +tests (35 passed, three intentional ignored), Clippy all targets/features, +CPU-only ICP, optional ICP WASM and formatting checks passed. + +## Completed-call performance + +Two separate sequential native release surveys completed successfully on +Vulkan with a discrete device, without recording the adapter model. Each +survey covers N=257/4097/65537/262145 source points and N+17 target points. +U means unique distances, G partially gated points, and D duplicate/tied +points; their trim fractions are 0.5/0.8/0.9 respectively. Every alignment +performs exactly three updates with zero stopping epsilons and a 0.25 gate. +The same input poses and cloud generation are used for all routes. + +Each of six stages receives three completed warmups and 20 paired rounds. +The stage order rotates each round: CPU; legacy GPU correspondences with +CPU trimming/moments; explicit CompensatedF32 resident reduction; default +DoubleSingle resident reduction; prepared F32; prepared DoubleSingle. +Compilation, input generation and independent validation occur outside the +timers. No local builds, tests or statistical analysis ran concurrently +with either survey. Both runs use already warmed pipeline caches; cold +shader compilation and construction of the runtime are excluded. + +Complete alignment timers include CPU packing/tree construction, explicit +uploads, GPU work/completion, fixed-size readbacks, CPU Horn solve and +composition, history and final all-gated fitness. Prepared timers use one +already uploaded input and the SAME final DoubleSingle pose for both +precisions; they include submission and wait_idle completion, and exclude +input preparation/upload, output readback and the CPU pose solve. Prepared +outputs are downloaded and checked after timing. Every timed call also +matches its own full warm control result exactly. Counts, independent f64 +fitness and the unchanged 3e-5 precise CPU pose threshold are checked for +all profiles, including the previously adverse large gated case. + +Times below are medians in milliseconds. Ratios are the geometric mean of +20 paired baseline/DoubleSingle time ratios; values above one mean the +precise route is faster than that baseline. Parenthesized intervals use +10,000 paired log-ratio bootstrap draws (indices 249/9749), seeded with +20261013 + 3*record_index + baseline_index. They describe within-survey +variation only. Scheduling changed absolute times between surveys; these +synthetic, single-backend measurements do not establish universal speedups. + +| Survey | Source/target | Profile/trim | CPU ms | Legacy GPU ms | F32 ms | DoubleSingle ms | Prepared F32/DS ms | F32/DS pose error | CPU/DS ratio (95%) | Legacy/DS ratio (95%) | F32/DS ratio (95%) | +| --- | ---: | :---: | ---: | ---: | ---: | ---: | ---: | ---: | --- | --- | --- | +| 1 | 257/274 | U/0.5 | 0.285 | 7.998 | 5.553 | 5.384 | 0.702/0.591 | 5.96e-08/0 | 0.05x (0.05-0.06) | 1.49x (1.40-1.60) | 1.00x (0.94-1.06) | +| 1 | 257/274 | G/0.8 | 0.320 | 9.696 | 5.769 | 5.827 | 0.703/0.642 | 2.38e-07/0 | 0.06x (0.05-0.06) | 1.49x (1.37-1.61) | 1.00x (0.94-1.07) | +| 1 | 257/274 | D/0.9 | 0.380 | 8.523 | 5.474 | 5.411 | 0.681/0.612 | 1.19e-07/0 | 0.07x (0.06-0.07) | 1.50x (1.41-1.60) | 0.98x (0.92-1.05) | +| 1 | 4,097/4,114 | U/0.5 | 7.293 | 10.957 | 6.783 | 6.848 | 0.838/1.650 | 0/0 | 0.97x (0.90-1.05) | 1.52x (1.43-1.61) | 0.96x (0.91-1.03) | +| 1 | 4,097/4,114 | G/0.8 | 9.075 | 11.479 | 6.968 | 7.236 | 0.853/0.964 | 9.54e-07/2.38e-07 | 1.21x (1.11-1.32) | 1.55x (1.46-1.64) | 0.96x (0.90-1.03) | +| 1 | 4,097/4,114 | D/0.9 | 9.395 | 11.748 | 7.728 | 7.839 | 0.913/0.964 | 9.54e-07/9.54e-07 | 1.23x (1.13-1.34) | 1.48x (1.39-1.58) | 0.95x (0.90-1.01) | +| 1 | 65,537/65,554 | U/0.5 | 166.437 | 65.097 | 29.881 | 36.442 | 3.127/4.915 | 1.98e-06/0 | 4.55x (4.31-4.85) | 1.79x (1.69-1.93) | 0.80x (0.75-0.86) | +| 1 | 65,537/65,554 | G/0.8 | 174.936 | 63.685 | 32.264 | 38.300 | 3.928/5.179 | 5.25e-06/4.66e-10 | 4.87x (4.57-5.20) | 1.68x (1.59-1.78) | 0.85x (0.80-0.89) | +| 1 | 65,537/65,554 | D/0.9 | 204.321 | 66.742 | 31.049 | 37.248 | 3.127/4.895 | 9.54e-07/1.43e-06 | 5.67x (5.28-6.09) | 1.86x (1.77-1.97) | 0.84x (0.81-0.88) | +| 1 | 262,145/262,162 | U/0.5 | 832.833 | 270.507 | 149.388 | 154.841 | 11.603/16.674 | 7.15e-07/0 | 6.13x (5.14-7.38) | 1.82x (1.58-2.08) | 0.97x (0.89-1.07) | +| 1 | 262,145/262,162 | G/0.8 | 1081.770 | 290.580 | 138.351 | 144.208 | 13.604/15.220 | 0.000815/7.63e-06 | 7.10x (6.25-8.23) | 2.00x (1.78-2.28) | 0.96x (0.84-1.11) | +| 1 | 262,145/262,162 | D/0.9 | 1036.076 | 275.925 | 123.198 | 136.274 | 12.489/14.292 | 5.72e-06/1.91e-06 | 7.59x (6.95-8.32) | 2.03x (1.86-2.24) | 0.90x (0.83-0.97) | +| 2 | 257/274 | U/0.5 | 0.622 | 25.688 | 15.749 | 16.893 | 1.098/1.261 | 5.96e-08/0 | 0.04x (0.03-0.04) | 1.59x (1.32-1.91) | 1.01x (0.92-1.11) | +| 2 | 257/274 | G/0.8 | 0.503 | 13.182 | 8.732 | 9.280 | 1.292/1.355 | 2.38e-07/0 | 0.05x (0.05-0.06) | 1.46x (1.39-1.53) | 0.98x (0.94-1.02) | +| 2 | 257/274 | D/0.9 | 0.566 | 12.493 | 8.976 | 9.110 | 1.440/1.446 | 1.19e-07/0 | 0.06x (0.06-0.06) | 1.39x (1.29-1.50) | 0.97x (0.91-1.04) | +| 2 | 4,097/4,114 | U/0.5 | 11.265 | 16.485 | 11.219 | 13.794 | 2.037/2.606 | 0/0 | 0.78x (0.70-0.86) | 1.21x (1.11-1.35) | 0.83x (0.79-0.87) | +| 2 | 4,097/4,114 | G/0.8 | 15.925 | 27.373 | 15.722 | 17.740 | 3.060/3.792 | 9.54e-07/2.38e-07 | 0.77x (0.72-0.83) | 1.37x (1.30-1.44) | 0.85x (0.80-0.91) | +| 2 | 4,097/4,114 | D/0.9 | 12.697 | 16.876 | 12.075 | 14.239 | 2.274/2.778 | 9.54e-07/9.54e-07 | 0.86x (0.77-0.94) | 1.20x (1.16-1.25) | 0.85x (0.83-0.89) | +| 2 | 65,537/65,554 | U/0.5 | 238.765 | 98.429 | 54.655 | 64.162 | 7.585/9.317 | 1.98e-06/0 | 3.80x (3.56-4.08) | 1.56x (1.48-1.65) | 0.86x (0.84-0.88) | +| 2 | 65,537/65,554 | G/0.8 | 264.383 | 83.826 | 38.227 | 40.462 | 2.684/3.642 | 5.25e-06/4.66e-10 | 5.87x (4.97-6.85) | 1.94x (1.74-2.16) | 0.90x (0.84-0.96) | +| 2 | 65,537/65,554 | D/0.9 | 261.298 | 88.903 | 59.278 | 67.621 | 8.030/8.679 | 9.54e-07/1.43e-06 | 4.28x (3.92-4.67) | 1.47x (1.34-1.61) | 0.88x (0.85-0.91) | +| 2 | 262,145/262,162 | U/0.5 | 604.312 | 208.490 | 103.102 | 109.633 | 11.129/12.643 | 7.15e-07/0 | 5.93x (5.51-6.40) | 2.03x (1.88-2.19) | 0.97x (0.91-1.04) | +| 2 | 262,145/262,162 | G/0.8 | 700.950 | 237.571 | 112.726 | 127.851 | 13.543/15.225 | 0.000815/7.63e-06 | 5.79x (5.17-6.43) | 1.86x (1.74-1.99) | 0.93x (0.89-0.97) | +| 2 | 262,145/262,162 | D/0.9 | 875.522 | 251.212 | 111.059 | 120.424 | 5.507/6.887 | 5.72e-06/1.91e-06 | 7.28x (6.51-8.14) | 2.06x (1.85-2.30) | 0.91x (0.81-1.02) | + +For 262,145 points, complete DoubleSingle calls are 5.79-7.59x faster than +CPU and 1.82-2.06x faster than legacy GPU correspondence calls across both +surveys. Precision has a cost: they take 1.03-1.11x the time of the explicit +previous resident F32 route by paired geometric ratios. For 65,537 points, +the CPU gain is 3.80-5.87x, with 1.12-1.25x the previous resident time. +For 4,097 points the CPU/precise ratio ranges from 0.77 to 1.23 and depends +on the profile/survey. At 257 points CPU is 14.7-27.5x faster. No automatic +CPU/GPU choice is introduced. The explicit F32 option leaves that arithmetic +tradeoff to the caller. + +The adverse gated 262,145-point pose difference is 0.0008154977 for the +previous arithmetic and 0.0000076293945 for DoubleSingle in both surveys. +The duplicate 65,537-point profile has a slightly larger pose difference +with DoubleSingle (0.0000014305115 versus 0.0000009536743); both meet the +same threshold. Improved moments do not imply a uniformly smaller pose +error for every geometry and selected-set trajectory. + +At N=262145/T=262162, three updates require four fixed-size readbacks: +legacy H2D/D2H are 34,604,528/14,680,120 payload bytes; both resident +precisions use 13,633,200/704. Precision adds no payload transfer. + +Private raw timing records and derived artifacts remain outside Git: + +- `C:/Users/rsasa/Workspace/SpatialRust/target/automation/gpu-icp-trimmed-precision-benchmark-1.log` +- `C:/Users/rsasa/Workspace/SpatialRust/target/automation/gpu-icp-trimmed-precision-benchmark-2.log` +- `C:/Users/rsasa/Workspace/SpatialRust/target/automation/gpu-icp-trimmed-precision-measurements.json` +- `C:/Users/rsasa/Workspace/SpatialRust/target/automation/analyze-gpu-icp-trimmed-precision.py` + +Reproduce from `C:/Users/rsasa/Workspace/SpatialRust`: compile with `--no-run` +first, then run this exact test alone twice, sequentially, with no overlapping +local builds or benchmarks: + +```text +cargo test -p spatialrust-registration --no-default-features --features register-icp-gpu --release --test gpu_icp_trimmed_precision_benchmark completed_trimmed_icp_precision_calls -- --include-ignored --exact --nocapture --test-threads=1 +``` + +Source references: + +- `C:/Users/rsasa/Workspace/SpatialRust/crates/spatialrust-gpu/src/kernels/icp_trim.rs` +- `C:/Users/rsasa/Workspace/SpatialRust/crates/spatialrust-gpu/src/shaders/icp_trim_moments.wgsl` +- `C:/Users/rsasa/Workspace/SpatialRust/crates/spatialrust-gpu/src/kernels/icp_trim_math_tests.rs` +- `C:/Users/rsasa/Workspace/SpatialRust/crates/spatialrust-gpu/src/kernels/icp_trim_precision_tests.rs` +- `C:/Users/rsasa/Workspace/SpatialRust/crates/spatialrust-gpu/tests/fixtures/icp_trimmed_f32_previous.wgsl` +- `C:/Users/rsasa/Workspace/SpatialRust/crates/spatialrust-registration/src/icp_reduction.rs` +- `C:/Users/rsasa/Workspace/SpatialRust/crates/spatialrust-registration/tests/gpu_icp_trimmed_precision.rs` +- `C:/Users/rsasa/Workspace/SpatialRust/crates/spatialrust-registration/tests/gpu_icp_trimmed_precision_benchmark.rs`