diff --git a/ext/examples/all_products.rs b/ext/examples/all_products.rs new file mode 100644 index 0000000000..17992ab1d0 --- /dev/null +++ b/ext/examples/all_products.rs @@ -0,0 +1,84 @@ +//! Benchmark helper: compute *all* left-multiplication-by-generator products on a resolved module. +//! +//! By default this extends every product map together via [`ExtAlgebra::extend_all_products`], which +//! batches the quasi-inverse solve at each output bidegree (one shared solve per bidegree). Set +//! `EXT_PER_MAP=1` to fall back to extending each map on its own instead, and compare timings. +//! +//! Extending a product map is the chain-map lift that repeatedly calls +//! [`ChainComplex::apply_quasi_inverse`], so this exercises the quasi-inverse path across the whole +//! plane. Set `EXT_NASSAU_NO_SAVE_QI=1` on the resolve to store only the differentials, or +//! `EXT_NASSAU_RECOMPUTE_QI=1` here to force recompute-on-demand even when quasi-inverses are saved. +//! +//! With `EXT_DUMP_PRODUCTS=1` it prints the full multiplication table (sorted) instead of a timing +//! line, so the batched and per-map paths can be diffed. + +use std::{sync::Arc, time::Instant}; + +use ext::{ + chain_complex::{ChainComplex, FreeChainComplex}, + ext_algebra::ExtAlgebra, + utils::query_module, +}; + +fn main() -> anyhow::Result<()> { + ext::utils::init_logging()?; + + // Loading (from the save dir) is untimed; only the product-map extension is timed. + let resolution = Arc::new(query_module(None, true)?); + let e2 = ExtAlgebra::from_resolution(resolution)?; + + let per_map = std::env::var_os("EXT_PER_MAP").is_some(); + + let start = Instant::now(); + let mut num_maps = 0usize; + if per_map { + // Old path: extend each product map on its own (recomputes each qi once per map). + for b in e2.resolution().iter_stem() { + for g in e2.basis(b) { + e2.generator_product_map(g).extend_all(); + num_maps += 1; + } + } + } else { + // Batched path: extend all product maps together in bidegree-major order (one shared qi + // solve per output bidegree). + num_maps = e2 + .resolution() + .iter_stem() + .map(|b| e2.resolution().number_of_gens_in_bidegree(b)) + .sum(); + e2.extend_all_products(); + } + let elapsed = start.elapsed(); + + let mode = if per_map { "per-map" } else { "batched" }; + eprintln!("all_products ({mode}): extended {num_maps} product maps in {elapsed:.3?}"); + + // Correctness dump: the full multiplication table, sorted so batched/per-map can be diffed. + if std::env::var_os("EXT_DUMP_PRODUCTS").is_some() { + let mut lines = Vec::new(); + for mb in e2.resolution().iter_stem() { + for mg in e2.basis(mb) { + let x = e2.generator(mg); + for b in e2.unit().iter_nonzero_stem() { + let Some(rows) = e2.multiply_into(&x, b) else { + continue; + }; + for (g, row) in e2.unit_basis(b).into_iter().zip(rows.iter()) { + let coords: Vec = row.iter().collect(); + if coords.iter().any(|&c| c != 0) { + lines.push(format!("x_{mg} · x_{g} = {coords:?}")); + } + } + } + } + } + lines.sort(); + for l in lines { + println!("{l}"); + } + } else { + println!("{num_maps} {}", elapsed.as_secs_f64()); + } + Ok(()) +} diff --git a/ext/src/chain_complex/chain_homotopy.rs b/ext/src/chain_complex/chain_homotopy.rs index c0ec0fd0d2..7d77d95882 100644 --- a/ext/src/chain_complex/chain_homotopy.rs +++ b/ext/src/chain_complex/chain_homotopy.rs @@ -11,7 +11,7 @@ use sseq::coordinates::{Bidegree, BidegreeRange}; use crate::{ chain_complex::{ChainComplex, FreeChainComplex}, - resolution_homomorphism::ResolutionHomomorphism, + resolution_homomorphism::{LiftRequest, Liftable, ResolutionHomomorphism}, save::{SaveDirectory, SaveKind}, }; @@ -153,12 +153,37 @@ impl< } fn extend_step(&self, source: Bidegree) -> std::ops::Range { + match self.prepare_step(source) { + HomotopyPrep::Done(range) => range, + HomotopyPrep::NeedsLift(mut pending) => { + let p = self.prime(); + let mut outputs = + vec![FpVector::new(p, pending.target_dim); pending.scratches.len()]; + assert!(U::apply_quasi_inverse( + &*self.right.target, + &mut outputs, + pending.target, + &pending.scratches, + )); + self.finish_step(&mut pending, &outputs) + } + } + } + + /// The part of a homotopy step *before* the quasi-inverse solve. Resolves the cases that need + /// no quasi-inverse (already computed, zero, or loaded from the save store) itself, returning + /// [`HomotopyPrep::Done`]; otherwise returns [`HomotopyPrep::NeedsLift`] with the vectors to + /// lift at `target = source + (1, 0) - shift`, to be completed by [`Self::finish_step`]. + /// + /// Splitting the step this way lets [`MultiLift`](crate::resolution_homomorphism::MultiLift) + /// batch the quasi-inverse solve of many homotopies at a shared bidegree. + fn prepare_step(&self, source: Bidegree) -> HomotopyPrep { let p = self.prime(); let shift = self.shift(); let target = source + Bidegree::s_t(1, 0) - shift; if self.homotopies[source.s()].next_degree() > source.t() { - return source.t()..source.t() + 1; + return HomotopyPrep::Done(source.t()..source.t() + 1); } let num_gens = self @@ -175,7 +200,9 @@ impl< // these values. if target.s() == 0 || target_dim == 0 || num_gens == 0 { let outputs = vec![FpVector::new(p, target_dim); num_gens]; - return self.homotopies[source.s()].add_generators_from_rows_ooo(source.t(), outputs); + return HomotopyPrep::Done( + self.homotopies[source.s()].add_generators_from_rows_ooo(source.t(), outputs), + ); } if let Some(dir) = self.save_dir.read() @@ -189,11 +216,11 @@ impl< for _ in 0..num_gens { outputs.push(FpVector::from_bytes(p, target_dim, &mut f).unwrap()); } - return self.homotopies[source.s()].add_generators_from_rows_ooo(source.t(), outputs); + return HomotopyPrep::Done( + self.homotopies[source.s()].add_generators_from_rows_ooo(source.t(), outputs), + ); } - let mut outputs = vec![FpVector::new(p, target_dim); num_gens]; - let f = |i| { let mut scratch = FpVector::new( p, @@ -256,24 +283,32 @@ impl< let scratches: Vec = (0..num_gens).into_maybe_par_iter().map(f).collect(); - assert!(U::apply_quasi_inverse( - &*self.right.target, - &mut outputs, + HomotopyPrep::NeedsLift(PendingHomotopy { + source, target, - &scratches, - )); + target_dim, + scratches, + }) + } + /// Complete a homotopy step: persist and register the lifted `outputs`. + fn finish_step( + &self, + pending: &mut PendingHomotopy, + outputs: &[FpVector], + ) -> std::ops::Range { if let Some(dir) = self.save_dir.write() { let mut f = self .left .source - .save_file(SaveKind::ChainHomotopy, source) + .save_file(SaveKind::ChainHomotopy, pending.source) .create_file(dir.to_owned(), false); - for row in &outputs { + for row in outputs { row.to_bytes(&mut f).unwrap(); } } - self.homotopies[source.s()].add_generators_from_rows_ooo(source.t(), outputs) + self.homotopies[pending.source.s()] + .add_generators_from_rows_ooo(pending.source.t(), outputs.to_vec()) } pub fn homotopy(&self, source_s: i32) -> Arc> { @@ -285,6 +320,59 @@ impl< } } +/// Outcome of [`ChainHomotopy::prepare_step`]. +enum HomotopyPrep { + /// The step needed no quasi-inverse and is already finished; carries the range of + /// newly-contiguous source degrees. + Done(std::ops::Range), + /// The step needs a quasi-inverse solve at `target`; complete it with + /// [`ChainHomotopy::finish_step`]. + NeedsLift(PendingHomotopy), +} + +/// A homotopy step awaiting its quasi-inverse solve; see [`ChainHomotopy::prepare_step`]. +struct PendingHomotopy { + source: Bidegree, + target: Bidegree, + target_dim: usize, + /// The vectors to lift, one per source generator. + scratches: Vec, +} + +impl< + S: FreeChainComplex + Send + Sync + 'static, + T: FreeChainComplex + Send + Sync + 'static, + U: ChainComplex + Send + Sync + 'static, +> Liftable for ChainHomotopy +{ + fn prepare(&self, b: Bidegree) -> Option> { + // `b` is the target bidegree where the quasi-inverse is taken; recover the source step. + let shift = self.shift(); + let source = b + shift - Bidegree::s_t(1, 0); + // Below the bottom homotopy row, or out of the source/target computed range: no work here. + // (The zero row `target.s() == 0` is *not* skipped — `prepare_step` finishes it itself.) + if source.s() < shift.s() - 1 || !self.left.source.has_computed_bidegree(source) { + return None; + } + if b.s() >= 1 && !self.right.target.has_computed_bidegree(b) { + return None; + } + self.initialize_homotopies(source.s() + 1); + match self.prepare_step(source) { + HomotopyPrep::Done(_) => None, + HomotopyPrep::NeedsLift(mut pending) => { + let inputs = std::mem::take(&mut pending.scratches); + Some(LiftRequest { + inputs, + finish: Box::new(move |outputs| { + self.finish_step(&mut pending, outputs); + }), + }) + } + } + } +} + // The secondary lift of a `ChainHomotopy` lives here, beside the primary object it lifts, rather // than in the monolithic `secondary` module. This keeps `secondary.rs` to the shared lift machinery // and pairs each variant with its primary for locality. The module is `pub(crate)`; diff --git a/ext/src/ext_algebra/massey.rs b/ext/src/ext_algebra/massey.rs index d29acdc403..30508fd5d3 100644 --- a/ext/src/ext_algebra/massey.rs +++ b/ext/src/ext_algebra/massey.rs @@ -27,7 +27,7 @@ use sseq::coordinates::{Bidegree, BidegreeElement, BidegreeGenerator}; use super::ExtAlgebra; use crate::{ chain_complex::{AugmentedChainComplex, ChainHomotopy, FreeChainComplex}, - resolution_homomorphism::ResolutionHomomorphism, + resolution_homomorphism::{Liftable, MultiLift, ResolutionHomomorphism}, }; /// The result of a Massey product computation @@ -145,11 +145,6 @@ where let representative = if target_num_gens == 0 { FpVector::new(p, 0) } else { - // Where `a`'s generators sit in the homotopy output, so we can pair against them. - let offset_a = - unit.module(a.degree().s()) - .generator_offset(a.degree().t(), a.degree().t(), 0); - let a_coords: Vec = a.vec().iter().collect(); let c_coords: Vec = c.vec().iter().collect(); let f_c = Arc::new(ResolutionHomomorphism::from_class( @@ -164,17 +159,7 @@ where let homotopy = ChainHomotopy::new(f_c, b_hom); homotopy.extend(tot); - let last = homotopy.homotopy(tot.s()); - let mut representative = FpVector::new(p, target_num_gens); - for i in 0..target_num_gens { - let output = last.output(tot.t(), i); - for (k, &val) in a_coords.iter().enumerate() { - if val != 0 { - representative.add_basis_element(i, val * output.entry(offset_a + k)); - } - } - } - representative + self.massey_read_representative(a, &homotopy, tot) }; let indeterminacy = self.massey_indeterminacy(a, c, tot); @@ -184,6 +169,39 @@ where }) } + /// Read the bracket representative off an already-extended null-homotopy of `b ∘ c` (the + /// `homotopy`), by pairing its top level (filtration `tot.s()`) against the first factor `a`. + /// + /// Factored out of [`massey_bracket_of`](Self::massey_bracket_of) so the batched + /// [`massey_iter_c`](Self::massey_iter_c), which extends all the null-homotopies together through + /// [`MultiLift`] before reading any bracket, shares the exact same read. + fn massey_read_representative( + &self, + a: &BidegreeElement, + homotopy: &ChainHomotopy, + tot: Bidegree, + ) -> FpVector { + let p = self.prime(); + let target_num_gens = self.resolution().number_of_gens_in_bidegree(tot); + // Where `a`'s generators sit in the homotopy output, so we can pair against them. + let offset_a = + self.unit() + .module(a.degree().s()) + .generator_offset(a.degree().t(), a.degree().t(), 0); + let a_coords: Vec = a.vec().iter().collect(); + let last = homotopy.homotopy(tot.s()); + let mut representative = FpVector::new(p, target_num_gens); + for i in 0..target_num_gens { + let output = last.output(tot.t(), i); + for (k, &val) in a_coords.iter().enumerate() { + if val != 0 { + representative.add_basis_element(i, val * output.entry(offset_a + k)); + } + } + } + representative + } + /// Compute a representative of a Massey product evaluated at `row` from the per-generator /// bracket matrix `answers`. Used by [`massey_iter_a`](Self::massey_iter_a), which builds one /// null-homotopy for fixed `b, c` and reads a whole family of first factors off `answers`. @@ -245,32 +263,103 @@ where /// /// Brackets that contain `0` are omitted. /// - /// This iterates over the third factor, building a fresh null-homotopy of `b ∘ c` per `c`. To - /// vary the *first* factor with `b, c` fixed instead, use - /// [`massey_iter_a`](Self::massey_iter_a). + /// This iterates over the third factor. Rather than building and extending a fresh null-homotopy + /// of `b ∘ c` per `c` — which, with quasi-inverses recomputed on demand, would re-solve the + /// unit's quasi-inverse at every bidegree once per third factor — it builds every map up front + /// and extends them together through [`MultiLift`], so each quasi-inverse is solved once and + /// shared across all third factors. To vary the *first* factor with `b, c` fixed instead, use + /// [`massey_iter_a`](Self::massey_iter_a), which already reuses a single null-homotopy. pub fn massey_iter_c( &self, a: &BidegreeElement, b: &BidegreeElement, - ) -> Vec<(BidegreeElement, MasseyResult)> { + ) -> Vec<(BidegreeElement, MasseyResult)> + where + CC: 'static, + { + let resolution = self.resolution(); + let unit = self.unit(); let shift = Self::massey_shift(a, b); - let b_hom = self.massey_b_hom(b, shift); - let mut results = Vec::new(); - for c_deg in self.resolution().iter_nonzero_stem() { + // Multiplication-by-`b` self-map of the unit. Built unextended so it can ride along in the + // first batch below (it shares the unit target with every `f_c`). + let b_coords: Vec = b.vec().iter().collect(); + let b_hom = Arc::new(ResolutionHomomorphism::from_class( + String::new(), + Arc::clone(unit), + Arc::clone(unit), + b.degree(), + &b_coords, + )); + + // Enumerate every valid third factor `c` (kernel of `· b`) whose bracket lands in a computed, + // non-empty bidegree, and build its realising map `f_c` and null-homotopy of `b ∘ c`. + // Brackets landing in an empty group are the (defined) zero element — they always contain + // zero, so they are dropped here exactly as the final `contains_zero` filter would. + let mut pending: Vec<(BidegreeElement, Bidegree, Arc>)> = + Vec::new(); + let mut f_cs: Vec>> = Vec::new(); + for c_deg in resolution.iter_nonzero_stem() { let Some(kernel) = self.massey_kernel(b, c_deg) else { continue; }; + let tot = c_deg + shift; + if !resolution.has_computed_bidegree(tot) + || resolution.number_of_gens_in_bidegree(tot) == 0 + { + continue; + } for row in kernel.iter() { let c = BidegreeElement::new(c_deg, row.to_owned()); - let Some(result) = self.massey_bracket_of(a, Arc::clone(&b_hom), shift, &c) else { - continue; - }; - if result.contains_zero() { - continue; - } - results.push((c, result)); + let c_coords: Vec = c.vec().iter().collect(); + let f_c = Arc::new(ResolutionHomomorphism::from_class( + String::new(), + Arc::clone(resolution), + Arc::clone(unit), + c_deg, + &c_coords, + )); + let homotopy = Arc::new(ChainHomotopy::new(Arc::clone(&f_c), Arc::clone(&b_hom))); + f_cs.push(f_c); + pending.push((c, tot, homotopy)); + } + } + + // No valid third factors: nothing to read, and in particular no reason to extend `b_hom`. + if pending.is_empty() { + return Vec::new(); + } + + // First batch: `b_hom` and every `f_c` against the unit. The homotopies read these maps, so + // they must be fully extended before the homotopies are lifted. Every one of these maps is + // read only up to output bidegree `shift` (the map value at `tot = c_deg + shift` lands at + // `tot - c_deg = shift`), so bound the sweep there rather than across the whole plane. + let mut map_liftables: Vec> = Vec::with_capacity(f_cs.len() + 1); + map_liftables.push(Arc::clone(&b_hom) as Arc); + map_liftables.extend(f_cs.iter().map(|f| Arc::clone(f) as Arc)); + MultiLift::new(Arc::clone(unit), map_liftables).extend_through_stem(shift); + + // Second batch: every null-homotopy against the unit. Each bracket is read at `tot`, i.e. + // the homotopy's top lands at output bidegree `a.degree()`, so bound the sweep there. + let homotopy_liftables: Vec> = pending + .iter() + .map(|(_, _, h)| Arc::clone(h) as Arc) + .collect(); + MultiLift::new(Arc::clone(unit), homotopy_liftables).extend_through_stem(a.degree()); + + // Read off each bracket from its extended homotopy. + let mut results = Vec::new(); + for (c, tot, homotopy) in pending { + let representative = self.massey_read_representative(a, &homotopy, tot); + let indeterminacy = self.massey_indeterminacy(a, &c, tot); + let result = MasseyResult { + degree: tot, + coset: AffineSubspace::new(representative, indeterminacy), + }; + if result.contains_zero() { + continue; } + results.push((c, result)); } results } diff --git a/ext/src/ext_algebra/mod.rs b/ext/src/ext_algebra/mod.rs index f229886636..1b4ae2a3a2 100644 --- a/ext/src/ext_algebra/mod.rs +++ b/ext/src/ext_algebra/mod.rs @@ -31,7 +31,7 @@ use sseq::coordinates::{Bidegree, BidegreeElement, BidegreeGenerator}; pub use self::secondary::{SecondaryExtAlgebra, SecondaryProduct}; use crate::{ chain_complex::{AugmentedChainComplex, FreeChainComplex}, - resolution_homomorphism::ResolutionHomomorphism, + resolution_homomorphism::{Liftable, MultiLift, ResolutionHomomorphism}, utils::{QueryModuleResolution, get_unit}, }; @@ -196,6 +196,33 @@ where Arc::clone(self.products.entry(g).or_insert(hom).value()) } + /// Build the product map for every generator of $\Ext(M, k)$ in the computed range and extend + /// them all together, in bidegree-major order, via [`MultiLift`]. + /// + /// All the product maps share the unit resolution as target, so lifting them together means the + /// unit's quasi-inverse at each output bidegree is solved once and reused across every product, + /// rather than once per product. With quasi-inverses recomputed on demand this keeps the cost of + /// computing all products at ~1x the on-disk-quasi-inverse baseline instead of scaling with the + /// number of generators. + pub fn extend_all_products(&self) + where + CC: Sync + 'static, + { + let liftables: Vec> = self + .resolution + .iter_stem() + .flat_map(|b| { + (0..self.resolution.number_of_gens_in_bidegree(b)) + .map(move |i| BidegreeGenerator::new(b, i)) + }) + .map(|g| self.generator_product_map(g) as Arc) + .collect(); + if liftables.is_empty() { + return; + } + MultiLift::new(Arc::clone(&self.unit), liftables).extend_all(); + } + /// Left-multiplication by the class `x` (in $\Ext(M, k)$), applied to every basis generator of /// $\Ext(k, k)$ at bidegree `b`. /// diff --git a/ext/src/ext_algebra/secondary.rs b/ext/src/ext_algebra/secondary.rs index 983a7c60ac..9c97998a62 100644 --- a/ext/src/ext_algebra/secondary.rs +++ b/ext/src/ext_algebra/secondary.rs @@ -17,7 +17,7 @@ use std::sync::{Arc, Mutex}; use algebra::pair_algebra::PairAlgebra; use dashmap::DashMap; use fp::{matrix::Subquotient, prime::Prime, vector::FpVector}; -use sseq::coordinates::{Bidegree, BidegreeElement}; +use sseq::coordinates::{Bidegree, BidegreeElement, BidegreeGenerator}; use super::ExtAlgebra; use crate::{ @@ -25,6 +25,7 @@ use crate::{ resolution_homomorphism::ResolutionHomomorphism, secondary::{ LAMBDA_BIDEGREE, SecondaryLift, SecondaryResolution, SecondaryResolutionHomomorphism, + batch_extend_secondary, }, }; @@ -208,6 +209,48 @@ where ) } + /// Build the secondary product lift for every generator of $\Ext(M, k)$ in the computed range + /// and extend them all together, in bidegree-major order, via [`MultiLift`](crate::resolution_homomorphism::MultiLift) — the secondary + /// analogue of [`ExtAlgebra::extend_all_products`](ExtAlgebra::extend_all_products). Must be + /// called after [`extend_all`](Self::extend_all). + /// + /// Every lift shares the unit resolution as target — for both the underlying primary maps and + /// the secondary homotopies — so lifting them together means the unit's quasi-inverse at each + /// bidegree is solved once and reused across every product rather than once per product. With + /// quasi-inverses recomputed on demand this keeps the cost of all the secondary products at ~1x + /// the on-disk baseline instead of scaling with the number of generators. + pub fn extend_all_secondary_products(&self) + where + CC: Send + Sync + 'static, + { + let resolution = self.alg.resolution(); + let lifts: Vec>> = resolution + .iter_stem() + // Skip filtration-zero classes (e.g. the unit): the secondary product lift there has + // `shift.s() == 1`, below where the source's secondary homotopies begin, so it is not a + // valid input to the secondary machinery (multiplying by such a class is degenerate). + .filter(|b| b.s() >= 1) + .flat_map(|b| { + (0..resolution.number_of_gens_in_bidegree(b)) + .map(move |i| BidegreeGenerator::new(b, i)) + }) + .map(|g| self.secondary_product_lift(&self.alg.generator(g))) + .collect(); + if lifts.is_empty() { + return; + } + + // Extend the underlying primary maps first (the secondary composites read them), each to its + // own extent exactly as `secondary_multiply_into` does — the secondary machinery's `max()` + // assumes that bound, so they are *not* over-extended through `MultiLift` here. Then batch + // the secondary homotopy solves, which all share the unit's quasi-inverse. + for lift in &lifts { + lift.underlying().extend_all(); + } + + batch_extend_secondary(Arc::clone(self.alg.unit()), lifts); + } + /// The secondary product of `x` with every $E_3$-surviving class of the unit at bidegree `b`, /// computed in $\Mod_{C\lambda^2}$. /// @@ -269,7 +312,7 @@ mod tests { use sseq::coordinates::BidegreeGenerator; use super::*; - use crate::utils::construct_standard; + use crate::{chain_complex::ChainComplex, utils::construct_standard}; #[test] fn test_sphere_d2() { @@ -305,4 +348,65 @@ mod tests { let h4_survives = sec_e2.survives(&h4).expect("h4 should have a computed d2"); assert!(!h4_survives, "h4 should not survive d2"); } + + /// Batching every secondary product together through [`MultiLift`](crate::resolution_homomorphism::MultiLift) + /// ([`extend_all_secondary_products`](SecondaryExtAlgebra::extend_all_secondary_products)) must + /// give the same $\Mod_{C\lambda^2}$ products as extending each lift on its own (the per-call + /// path inside [`secondary_multiply_into`](SecondaryExtAlgebra::secondary_multiply_into)). This + /// pins the multi-lift secondary batching end to end. + #[test] + fn batched_secondary_products_match() { + let build = || { + let res = Arc::new(construct_standard::("S_2", None).unwrap()); + res.compute_through_stem(Bidegree::n_s(14, 6)); + let e2 = Arc::new(ExtAlgebra::new(Arc::clone(&res), res)); + let sec = SecondaryExtAlgebra::new(Arc::clone(&e2)); + sec.extend_all(); + (e2, sec) + }; + let (e2_ref, sec_ref) = build(); + let (e2_bat, sec_bat) = build(); + + // Extend every secondary product together, in one batched pass. + sec_bat.extend_all_secondary_products(); + + let mut compared = 0; + for (n, s) in [(0, 1), (1, 1), (3, 1)] { + let g = BidegreeGenerator::new(Bidegree::n_s(n, s), 0); + let x_ref = e2_ref.generator(g); + let x_bat = e2_bat.generator(g); + + for b in e2_ref.resolution().iter_stem() { + // `secondary_multiply_into` reads up to `b + x.degree() + λ`, so only query where + // that lands in the computed range. + if b.n() > 8 + || !e2_ref + .resolution() + .has_computed_bidegree(b + x_ref.degree() + LAMBDA_BIDEGREE) + { + continue; + } + // Reference extends this lift on its own; batched reuses the pre-extended one. + let pr = sec_ref.secondary_multiply_into(&x_ref, b); + let pb = sec_bat.secondary_multiply_into(&x_bat, b); + assert_eq!( + pr.len(), + pb.len(), + "product count differs at x={x_ref}, b={b}" + ); + for (r, t) in pr.iter().zip(&pb) { + assert_eq!( + r.ext_part, t.ext_part, + "ext part differs at x={x_ref}, b={b}" + ); + assert_eq!( + r.lambda_part, t.lambda_part, + "lambda part differs at x={x_ref}, b={b}" + ); + compared += 1; + } + } + } + assert!(compared > 0, "expected to compare some secondary products"); + } } diff --git a/ext/src/nassau.rs b/ext/src/nassau.rs index 8d1919136a..20ecbf9858 100644 --- a/ext/src/nassau.rs +++ b/ext/src/nassau.rs @@ -15,7 +15,7 @@ use std::{ fmt::Display, io, - sync::{Arc, Mutex, mpsc}, + sync::{Arc, LazyLock, Mutex, mpsc}, }; use algebra::{ @@ -388,6 +388,17 @@ enum Magic { Fix = -3, } +/// Whether to persist quasi-inverses to disk during resolution. Disabled by +/// `EXT_NASSAU_NO_SAVE_QI`, in which case only the differentials are written (the quasi-inverses are +/// ~260-460x larger) and every downstream lift recomputes its quasi-inverse on demand. +static SAVE_QI: LazyLock = + LazyLock::new(|| std::env::var_os("EXT_NASSAU_NO_SAVE_QI").is_none()); + +/// Force `apply_quasi_inverse` to recompute rather than read a saved quasi-inverse, even when one +/// exists on disk. Used to measure the recompute cost in isolation. +static RECOMPUTE_QI: LazyLock = + LazyLock::new(|| std::env::var_os("EXT_NASSAU_RECOMPUTE_QI").is_some()); + /// A resolution of `S_2` using Nassau's algorithm. /// /// This aims to have an API similar to that of @@ -598,7 +609,9 @@ impl> Resolution { let next = &self.modules[b.s() - 2]; next.compute_basis(b.t()); - let mut f = if let Some(dir) = self.save_dir().write() { + // Skip writing the quasi-inverse when `EXT_NASSAU_NO_SAVE_QI` is set; `apply_quasi_inverse` + // recomputes it on demand from the differential (see `RecomputeReader`). + let mut f = if *SAVE_QI && let Some(dir) = self.save_dir().write() { let mut f = self .save_file(SaveKind::NassauQi, b - Bidegree::s_t(1, 0)) .create_file(dir.to_owned(), true); @@ -1057,14 +1070,23 @@ impl> ChainComplex for Resolution { for<'a> &'a mut T: Into>, for<'a> &'a S: Into>, { - let mut f = if let Some(dir) = self.save_dir.read() { - if let Some(f) = self.save_file(SaveKind::NassauQi, b).open_file(dir.clone()) { - f - } else { - return false; - } + // Read the saved quasi-inverse unless recomputation is forced. Fall back to recomputation + // whenever no saved stream is available: no store, the qis were never persisted + // (`EXT_NASSAU_NO_SAVE_QI`), or the store has no qi for this bidegree. The last case + // legitimately happens at the top of the computed region — nassau writes qi(s, t) while + // computing (s + 1, t), so qi(max_s, t) is never saved even though lifting into it is + // well-defined. `RecomputeReader` regenerates the exact same byte stream from + // `differentials[b.s]`, so the loop below is unchanged. + let saved = if *RECOMPUTE_QI { + None + } else if let Some(dir) = self.save_dir.read() { + self.save_file(SaveKind::NassauQi, b).open_file(dir.clone()) } else { - return false; + None + }; + let mut f: Box = match saved { + Some(f) => f, + None => Box::new(RecomputeReader::new(self, b)), }; let p = self.prime(); @@ -1224,6 +1246,148 @@ impl> ChainComplex for Resolution { } } +/// A streaming [`io::Read`] that regenerates the quasi-inverse of `d_{b.s}` at bidegree `b` on the +/// fly, in exactly the byte format that [`Resolution::write_qi`] wrote to disk. +/// +/// This lets [`ChainComplex::apply_quasi_inverse`] read a recomputed quasi-inverse through the same +/// code path as a saved one — the apply loop is unchanged; only the source of bytes differs. The +/// quasi-inverse is re-derived from `differentials[b.s]` alone (plus the module bases and the +/// deterministically-chosen subalgebra), never the rest of the resolution. +/// +/// It advances one signature at a time, holding only that signature's matrices, so its peak memory +/// matches what resolving this bidegree originally required — never the whole quasi-inverse, which +/// can reach hundreds of GB at record stems. +/// +/// Because the resolution is fully computed by the time a lift is requested, the recomputed +/// quasi-inverse always uses complete information, so it never emits a [`Magic::Fix`]. +struct RecomputeReader<'a, M: ZeroModule> { + res: &'a Resolution, + b: Bidegree, + subalgebra: MilnorSubalgebra, + algebra: Arc, + signatures: std::vec::IntoIter>, + scratch: FpVector, + buf: Vec, + pos: usize, + header_done: bool, + end_done: bool, +} + +impl<'a, M: ZeroModule> RecomputeReader<'a, M> { + fn new(res: &'a Resolution, b: Bidegree) -> Self { + let s = b.s(); + let t = b.t(); + // The subalgebra is chosen deterministically, as in `step_resolution_with_result` for the + // step that computed `(s + 1, t)` (the step that wrote this qi). + let subalgebra = MilnorSubalgebra::optimal_for( + Bidegree::s_t(s + 1, t) - Bidegree::s_t(0, res.max_degree), + ); + // `src` is the source of `d_s` (= F_s), `tgt` its target (= F_{s-1}). + let src = &res.modules[s]; + let tgt = &res.modules[s - 1]; + src.compute_basis(t); + tgt.compute_basis(t); + let algebra = tgt.algebra(); + + // Zero signature first, then the rest — matching the write order. + let signatures: Vec> = std::iter::once(subalgebra.zero_signature()) + .chain(subalgebra.iter_signatures(t)) + .collect(); + + Self { + res, + b, + subalgebra, + algebra, + signatures: signatures.into_iter(), + scratch: FpVector::new(res.prime(), 0), + buf: Vec::new(), + pos: 0, + header_done: false, + end_done: false, + } + } + + /// Write the qi header: target dimension, zero-signature masked source dimension, subalgebra. + fn write_header(&mut self) -> io::Result<()> { + let (s, t) = (self.b.s(), self.b.t()); + let target_dim = self.res.modules[s - 1].dimension(t); + let zero_mask_dim = self + .subalgebra + .signature_mask( + &self.algebra, + &self.res.modules[s], + t, + &self.subalgebra.zero_signature(), + ) + .count(); + self.buf.write_u64::(target_dim as u64)?; + self.buf.write_u64::(zero_mask_dim as u64)?; + self.subalgebra.to_bytes(&mut self.buf) + } + + /// Row-reduce `d_s` restricted to `signature` and append the resulting commands, exactly as + /// `write_qi` did during resolution. Writes nothing if the block has no pivots. + fn write_signature(&mut self, signature: &[PPartEntry]) -> io::Result<()> { + let (s, t) = (self.b.s(), self.b.t()); + let p = self.res.prime(); + + let src_mask: Vec = self + .subalgebra + .signature_mask(&self.algebra, &self.res.modules[s], t, signature) + .collect(); + let tgt_mask: Vec = self + .subalgebra + .signature_mask(&self.algebra, &self.res.modules[s - 1], t, signature) + .collect(); + + let full_matrix = { + let _guard = ParallelGuard::new(); + self.res.differentials[s].get_partial_matrix(t, &src_mask) + }; + let mut masked_matrix = + AugmentedMatrix::new(p, src_mask.len(), [tgt_mask.len(), src_mask.len()]); + masked_matrix + .segment(0, 0) + .add_masked(&full_matrix, &tgt_mask); + masked_matrix.segment(1, 1).add_identity(); + masked_matrix.row_reduce(); + + Resolution::::write_qi( + &mut Some(&mut self.buf), + &mut self.scratch, + signature, + &tgt_mask, + &full_matrix, + &masked_matrix, + ) + } +} + +impl> io::Read for RecomputeReader<'_, M> { + fn read(&mut self, out: &mut [u8]) -> io::Result { + while self.pos >= self.buf.len() { + self.buf.clear(); + self.pos = 0; + if !self.header_done { + self.write_header()?; + self.header_done = true; + } else if let Some(signature) = self.signatures.next() { + self.write_signature(&signature)?; + } else if !self.end_done { + self.buf.write_u64::(Magic::End as u64)?; + self.end_done = true; + } else { + return Ok(0); + } + } + let n = std::cmp::min(out.len(), self.buf.len() - self.pos); + out[..n].copy_from_slice(&self.buf[self.pos..self.pos + n]); + self.pos += n; + Ok(n) + } +} + impl> AugmentedChainComplex for Resolution { type ChainMap = FreeModuleHomomorphism; type TargetComplex = FiniteChainComplex>; diff --git a/ext/src/resolution_homomorphism.rs b/ext/src/resolution_homomorphism.rs index 271c7318d3..7d28b768dc 100644 --- a/ext/src/resolution_homomorphism.rs +++ b/ext/src/resolution_homomorphism.rs @@ -15,8 +15,8 @@ use fp::{ vector::{FpSliceMut, FpVector}, }; use maybe_rayon::prelude::*; -use once::OnceBiVec; -use sseq::coordinates::{Bidegree, BidegreeGenerator, BidegreeRange}; +use once::{OnceBiVec, OnceVec}; +use sseq::coordinates::{Bidegree, BidegreeGenerator, BidegreeRange, iter_s_t}; use crate::{ chain_complex::{AugmentedChainComplex, BoundedChainComplex, ChainComplex, FreeChainComplex}, @@ -201,12 +201,48 @@ where /// The user should call this function explicitly to manually define the chain map where the /// chain complex is not exact, and then call [`MuResolutionHomomorphism::extend_all`] to extend /// the rest by exactness. + /// + /// This is the single-map path: it recomputes/reads its own quasi-inverse inline. Callers + /// extending *several* maps with a common target should instead drive them through + /// [`MultiLift`], which batches the quasi-inverse solve at each output bidegree so it is + /// computed once and shared across all maps. #[tracing::instrument(skip(self, extra_images), fields(self = self.name, %input))] pub fn extend_step_raw( &self, input: Bidegree, extra_images: Option>, ) -> Range { + match self.prepare_step(input, extra_images) { + StepPrep::Done(range) => range, + StepPrep::NeedsLift(mut pending) => { + let p = self.source.prime(); + let mut results = + vec![FpVector::new(p, pending.fx_dimension); pending.fdx_vectors.len()]; + if !pending.fdx_vectors.is_empty() { + assert!(self.target.apply_quasi_inverse( + &mut results, + pending.output, + &pending.fdx_vectors + )); + } + self.finish_step(&mut pending, &results) + } + } + } + + /// The part of a lift step *before* the quasi-inverse solve: it resolves the cases that need no + /// quasi-inverse (already computed, loaded from the save store, zero-dimensional, or an + /// augmentation step) itself, returning [`StepPrep::Done`]; otherwise it returns + /// [`StepPrep::NeedsLift`] carrying the fdx vectors to lift at `output = input - shift` and the + /// partially-filled outputs, to be completed by [`Self::finish_step`]. + /// + /// Splitting the step this way lets [`MultiLift`] gather the fdx vectors of many maps at a + /// common output bidegree and issue a single batched `apply_quasi_inverse`. + pub(crate) fn prepare_step( + &self, + input: Bidegree, + extra_images: Option>, + ) -> StepPrep { let output = input - self.shift; assert!(self.target.has_computed_bidegree(output)); assert!(self.source.has_computed_bidegree(input)); @@ -216,7 +252,7 @@ where if input.t() < f_cur.next_degree() { assert!(extra_images.is_none()); // We need to signal to compute the dependents of this - return input.t()..input.t() + 1; + return StepPrep::Done(input.t()..input.t() + 1); } let p = self.source.prime(); @@ -225,45 +261,33 @@ where let fx_dimension = f_cur.target().dimension(output.t()); if num_gens == 0 || fx_dimension == 0 { - return f_cur.add_generators_from_rows_ooo( + return StepPrep::Done(f_cur.add_generators_from_rows_ooo( input.t(), vec![FpVector::new(p, fx_dimension); num_gens], - ); + )); } - if let Some(dir) = self.save_dir.read() { - let mut outputs = Vec::with_capacity(num_gens); - - if let Some(mut f) = self + if let Some(dir) = self.save_dir.read() + && let Some(mut f) = self .source .save_file(SaveKind::ChainMap, input) .open_file(dir.to_owned()) - { - let fx_dimension = f.read_u64::().unwrap() as usize; - for _ in 0..num_gens { - outputs.push(FpVector::from_bytes(p, fx_dimension, &mut f).unwrap()); - } - return f_cur.add_generators_from_rows_ooo(input.t(), outputs); + { + let fx_dimension = f.read_u64::().unwrap() as usize; + let mut outputs = Vec::with_capacity(num_gens); + for _ in 0..num_gens { + outputs.push(FpVector::from_bytes(p, fx_dimension, &mut f).unwrap()); } + return StepPrep::Done(f_cur.add_generators_from_rows_ooo(input.t(), outputs)); } if output.s() == 0 { let outputs = extra_images.unwrap_or_else(|| vec![FpVector::new(p, fx_dimension); num_gens]); - - if let Some(dir) = self.save_dir.write() { - let mut f = self - .source - .save_file(SaveKind::ChainMap, input) - .create_file(dir.clone(), false); - f.write_u64::(fx_dimension as u64).unwrap(); - for row in &outputs { - row.to_bytes(&mut f).unwrap(); - } - } - - return f_cur.add_generators_from_rows_ooo(input.t(), outputs); + self.save_chain_map(input, fx_dimension, &outputs); + return StepPrep::Done(f_cur.add_generators_from_rows_ooo(input.t(), outputs)); } + let mut outputs = vec![FpVector::new(p, fx_dimension); num_gens]; let d_source = self.source.differential(input.s()); let d_target = self.target.differential(output.s()); @@ -284,10 +308,14 @@ where } } - // Now do the rest + // Now compute the fdx vectors to lift, and remember which output row each fills. d_target.compute_auxiliary_data_through_degree(output.t()); - let compute_fdx_vector = |k| { + let qi_rows: Vec = (0..num_gens) + .filter(|&k| !d_source.output(input.t(), k).is_zero()) + .collect(); + + let compute_fdx_vector = |k: usize| { let dx_vector = d_source.output(input.t(), k); if dx_vector.is_zero() { None @@ -303,41 +331,240 @@ where } }; + // Same order as `qi_rows` (ascending generator index over non-zero differentials). let fdx_vectors: Vec = (0..num_gens) .into_maybe_par_iter() .filter_map(compute_fdx_vector) .collect(); - let mut qi_outputs: Vec<_> = outputs - .iter_mut() - .enumerate() - .filter_map(|(k, v)| { - if d_source.output(input.t(), k).is_zero() { - None - } else { - Some(v.as_slice_mut()) - } - }) - .collect(); + StepPrep::NeedsLift(PendingStep { + input, + output, + outputs, + fdx_vectors, + qi_rows, + fx_dimension, + }) + } - if !fdx_vectors.is_empty() { - assert!( - self.target - .apply_quasi_inverse(&mut qi_outputs, output, &fdx_vectors) - ); + /// Complete a lift step: scatter the lifted `results` into the pending outputs, persist the + /// chain map, and register the generators. `results[i]` is the lift of `pending.fdx_vectors[i]` + /// and fills output row `pending.qi_rows[i]`. + pub(crate) fn finish_step( + &self, + pending: &mut PendingStep, + results: &[FpVector], + ) -> Range { + assert_eq!(results.len(), pending.qi_rows.len()); + for (&k, result) in pending.qi_rows.iter().zip(results) { + pending.outputs[k].assign(result); } + self.save_chain_map(pending.input, pending.fx_dimension, &pending.outputs); + self.get_map(pending.input.s()) + .add_generators_from_rows_ooo(pending.input.t(), std::mem::take(&mut pending.outputs)) + } + fn save_chain_map(&self, input: Bidegree, fx_dimension: usize, outputs: &[FpVector]) { if let Some(dir) = self.save_dir.write() { let mut f = self .source .save_file(SaveKind::ChainMap, input) .create_file(dir.clone(), false); f.write_u64::(fx_dimension as u64).unwrap(); - for row in &outputs { + for row in outputs { row.to_bytes(&mut f).unwrap(); } } - f_cur.add_generators_from_rows_ooo(input.t(), outputs) + } +} + +/// Outcome of [`MuResolutionHomomorphism::prepare_step`]. +pub(crate) enum StepPrep { + /// The step needed no quasi-inverse and is already finished; carries the range of + /// newly-contiguous input degrees (as [`MuResolutionHomomorphism::extend_step_raw`] returns). + Done(Range), + /// The step needs a quasi-inverse solve at `output`; complete it with + /// [`MuResolutionHomomorphism::finish_step`]. + NeedsLift(PendingStep), +} + +/// A lift step awaiting its quasi-inverse solve; see [`MuResolutionHomomorphism::prepare_step`]. +pub(crate) struct PendingStep { + input: Bidegree, + output: Bidegree, + /// Outputs with the extra-image rows filled and the quasi-inverse rows left zero. + outputs: Vec, + /// The vectors to lift (`d(f(x))` preimages), one per non-zero-differential generator. + fdx_vectors: Vec, + /// `qi_rows[i]` is the output row that `fdx_vectors[i]`'s lift fills. + qi_rows: Vec, + /// Dimension of each lifted result (the target module dimension at `output`). + fx_dimension: usize, +} + +/// Something built by a sequence of quasi-inverse solves against a common target complex, one +/// target bidegree at a time. +/// +/// Implemented by chain-map extension ([`MuResolutionHomomorphism`]), the chain-homotopy lifts of +/// the Massey machinery ([`ChainHomotopy`](crate::chain_complex::ChainHomotopy)), and the secondary +/// lifts (via [`batch_extend_secondary`](crate::secondary::batch_extend_secondary)). They all lift +/// through the *same* target quasi-inverse at a given bidegree, so [`MultiLift`] can gather a batch +/// across implementors of different kinds and solve it once. +/// +/// The interface is deliberately free of the target/source type parameters: `prepare` returns plain +/// vectors to lift plus a boxed continuation, so the driver never needs to name a liftable's +/// internal state. +pub trait Liftable: Sync + Send { + /// Prepare the lift at target bidegree `b`. Returns `None` if this liftable has no quasi-inverse + /// work at `b` — out of range, already computed, or a step (augmentation, zero-dimensional) it + /// finished itself. Otherwise returns the vectors to lift at `b` and a continuation that + /// finishes the step once their lifts are known. + fn prepare(&self, b: Bidegree) -> Option>; +} + +/// The inputs to lift at one bidegree together with a continuation to finish the step; see +/// [`Liftable::prepare`]. +pub struct LiftRequest<'a> { + /// Vectors to lift, i.e. the `inputs` passed to `apply_quasi_inverse` at this bidegree. + pub inputs: Vec, + /// Called exactly once with the lifted `results` (`results[i]` lifts `inputs[i]`) to complete + /// the step (scatter the results, record generators). + pub finish: Box, +} + +impl Liftable for MuResolutionHomomorphism +where + CC1: FreeChainComplex, + CC1::Algebra: MuAlgebra, + CC2: ChainComplex, +{ + fn prepare(&self, b: Bidegree) -> Option> { + let input = b + self.shift; + if input.s() < self.shift.s() + || !self.source.has_computed_bidegree(input) + || !self.target.has_computed_bidegree(b) + { + return None; + } + match self.prepare_step(input, None) { + StepPrep::Done(_) => None, + StepPrep::NeedsLift(mut pending) => { + let inputs = std::mem::take(&mut pending.fdx_vectors); + Some(LiftRequest { + inputs, + finish: Box::new(move |results| { + self.finish_step(&mut pending, results); + }), + }) + } + } + } +} + +/// Extends several [`Liftable`]s that share one target complex, together, in bidegree-major order. +/// +/// At each output bidegree the inputs of every participating liftable are gathered and lifted with +/// a single [`ChainComplex::apply_quasi_inverse`], so the target's quasi-inverse there is solved +/// once and shared across all of them rather than recomputed once per liftable. This is what makes +/// recompute-on-demand (no saved quasi-inverses) cost ~1x across a many-map computation instead of +/// scaling with the number of maps. Concurrency across the plane is provided by +/// [`iter_s_t`], exactly as for a single map. +/// +/// Single-map callers do not need this — one map hits each bidegree once, so it already recomputes +/// each quasi-inverse once. Use [`MuResolutionHomomorphism::extend_all`] for those. +pub struct MultiLift { + target: Arc, + liftables: Vec>, +} + +impl MultiLift { + /// Build a driver over `liftables`, all of which must lift through `target`. + pub fn new(target: Arc, liftables: Vec>) -> Self { + Self { target, liftables } + } + + /// Extend every liftable as far as the shared target is resolved, batching the quasi-inverse + /// solve at each output bidegree. + pub fn extend_all(&self) { + self.extend_bounded(None); + } + + /// Like [`extend_all`](Self::extend_all), but only through the stem profile of `bound` — output + /// bidegrees `(s, t)` with `s <= bound.s()` and `n <= bound.n()`, the same shape + /// [`MuResolutionHomomorphism::extend_through_stem`] uses for a single map. Use this when the + /// batch's results are read only up to a known bidegree, so each liftable is extended to just + /// what it needs instead of across the whole computed plane. + pub fn extend_through_stem(&self, bound: Bidegree) { + self.extend_bounded(Some(bound)); + } + + /// Shared driver for [`extend_all`](Self::extend_all) and + /// [`extend_through_stem`](Self::extend_through_stem). `bound`, when present, caps the swept + /// output bidegrees to its stem profile (intersected with the target's computed range). + fn extend_bounded(&self, bound: Option) { + let mut max_s = self.target.next_homological_degree(); + if let Some(bound) = bound { + max_s = std::cmp::min(max_s, bound.s() + 1); + } + if max_s <= 0 || self.liftables.is_empty() { + return; + } + let min_t = self.target.min_degree(); + let min = Bidegree::s_t(0, min_t); + + // Per-output-row completion frontier, so `iter_s_t` can tell how far each row is done. Cell + // (s, t) is recorded at index `t - min_t` of row `s`; `push_ooo` returns the contiguous + // frontier that `iter_s_t` expects. + let completion: OnceVec> = OnceVec::new(); + for _ in 0..max_s { + completion.push(OnceVec::new()); + } + + let max_t = move |selff: &Self, s: i32| { + let mut t = selff.target.module(s).max_computed_degree() + 1; + if let Some(bound) = bound { + // Stem profile `n <= bound.n()`, i.e. `t <= bound.n() + s`, exclusive upper bound. + t = std::cmp::min(t, bound.n() + s + 1); + } + t + }; + let max = BidegreeRange::new(self, max_s, &max_t); + + iter_s_t(&|b| self.step_cell(b, &completion, min_t), min, max); + } + + /// Process one output bidegree: gather every participating liftable's inputs, do one batched + /// lift, then finish each. Returns the newly-contiguous frontier of this output row. + fn step_cell(&self, b: Bidegree, completion: &OnceVec>, min_t: i32) -> Range { + let p = self.target.prime(); + + let mut inputs: Vec = Vec::new(); + let mut finishers: Vec<(Box, Range)> = Vec::new(); + + for liftable in &self.liftables { + if let Some(mut req) = liftable.prepare(b) { + let start = inputs.len(); + inputs.append(&mut req.inputs); + finishers.push((req.finish, start..inputs.len())); + } + } + + // Every liftable at output `b` lifts through the same quasi-inverse, so all results have the + // target's module dimension at `b`. + let fx_dim = self.target.module(b.s()).dimension(b.t()); + let mut results = vec![FpVector::new(p, fx_dim); inputs.len()]; + if !inputs.is_empty() { + assert!(self.target.apply_quasi_inverse(&mut results, b, &inputs)); + } + // Run every finisher, even when a liftable contributed no inputs (a zero-dimensional step): + // its finish still has to register/extend the step so later reads of that bidegree see it. + // A no-input finisher receives an empty `results` slice. + for (finish, range) in finishers { + finish(&results[range]); + } + + let frontier = completion[b.s() as usize].push_ooo((), (b.t() - min_t) as usize); + (frontier.start as i32 + min_t)..(frontier.end as i32 + min_t) } } diff --git a/ext/src/secondary.rs b/ext/src/secondary.rs index e04b852a36..e5d0043d6d 100644 --- a/ext/src/secondary.rs +++ b/ext/src/secondary.rs @@ -34,6 +34,7 @@ pub use crate::{ }; use crate::{ chain_complex::{ChainComplex, FreeChainComplex}, + resolution_homomorphism::{LiftRequest, Liftable, MultiLift}, save::{SaveDirectory, SaveFile, SaveKind}, }; @@ -536,9 +537,39 @@ pub trait SecondaryLift: Sync + Sized { /// `self.try_compute_homotopy_step(b).unwrap()`. #[tracing::instrument(skip(self), fields(%b))] fn try_compute_homotopy_step(&self, b: Bidegree) -> anyhow::Result> { + match self.prepare_homotopy_step(b) { + SecondaryHomotopyPrep::Done(range) => Ok(range), + SecondaryHomotopyPrep::NeedsLift(pending) => { + let p = self.prime(); + let mut results = vec![FpVector::new(p, pending.target_dim); pending.num_gens]; + anyhow::ensure!( + self.target().apply_quasi_inverse( + &mut results, + pending.target_b, + &pending.intermediates, + ), + "secondary: failed to apply quasi-inverse at {b}; the input likely does not \ + lift" + ); + self.finish_homotopy_step(pending, &results) + } + } + } + + /// The part of a homotopy step *before* the quasi-inverse solve at `target_b = b - shift - + /// (0, 1)`. Resolves the cases that need no quasi-inverse (already computed, or loaded from the + /// save store) itself, returning [`SecondaryHomotopyPrep::Done`]; otherwise it assembles the + /// intermediates to lift and returns [`SecondaryHomotopyPrep::NeedsLift`], to be completed by + /// [`Self::finish_homotopy_step`]. + /// + /// Splitting the step this way lets [`MultiLift`] + /// gather the intermediates of many secondary lifts sharing a target at a common bidegree and + /// solve them with a single batched `apply_quasi_inverse`. Assumes the intermediates have + /// already been computed (via [`compute_intermediates`](Self::compute_intermediates)). + fn prepare_homotopy_step(&self, b: Bidegree) -> SecondaryHomotopyPrep { let homotopy = &self.homotopies()[b.s()]; if homotopy.homotopies.next_degree() > b.t() { - return Ok(b.t()..b.t() + 1); + return SecondaryHomotopyPrep::Done(b.t()..b.t() + 1); } let p = self.prime(); let shift = self.shift(); @@ -563,9 +594,11 @@ pub trait SecondaryLift: Sync + Sized { for _ in 0..num_gens { results.push(FpVector::from_bytes(p, target_dim, &mut f).unwrap()); } - return Ok(self.homotopies()[b.s()] - .homotopies - .add_generators_from_rows_ooo(b.t(), results)); + return SecondaryHomotopyPrep::Done( + self.homotopies()[b.s()] + .homotopies + .add_generators_from_rows_ooo(b.t(), results), + ); } } @@ -586,22 +619,52 @@ pub trait SecondaryLift: Sync + Sized { v }; - let mut intermediates: Vec = (0..num_gens) + let intermediates: Vec = (0..num_gens) .into_maybe_par_iter() .map(get_intermediate) .collect(); - let mut results = vec![FpVector::new(p, target_dim); num_gens]; + // The post-quasi-inverse lift-validity check (only at the bottom non-trivial row) consumes + // the *pre*-solve intermediates, but the batched driver moves them out to lift. Keep a copy + // for the check there; elsewhere it never runs. + let check = if b.s() == shift.s() + 1 { + Some(intermediates.clone()) + } else { + None + }; - anyhow::ensure!( - target.apply_quasi_inverse(&mut results, target_b, &intermediates,), - "secondary: failed to apply quasi-inverse at {b}; the input likely does not lift" - ); + SecondaryHomotopyPrep::NeedsLift(SecondaryPending { + b, + target_b, + num_gens, + target_dim, + intermediates, + check, + }) + } + + /// Complete a homotopy step: verify the lift (at the bottom non-trivial row), persist the + /// homotopy, delete the now-consumed intermediate files, and register the generators. + /// `results[i]` is the lift of the `i`th intermediate. See [`Self::prepare_homotopy_step`]. + fn finish_homotopy_step( + &self, + pending: SecondaryPending, + results: &[FpVector], + ) -> anyhow::Result> { + let SecondaryPending { + b, + target_b, + num_gens, + check, + .. + } = pending; + let p = self.prime(); - if b.s() == shift.s() + 1 { - // Check that we indeed had a lift - let d = target.differential(target_b.s()); - for (src, tgt) in std::iter::zip(&results, &mut intermediates) { + if b.s() == self.shift().s() + 1 { + // Check that we indeed had a lift. + let mut check = check.expect("check copy is populated at the bottom non-trivial row"); + let d = self.target().differential(target_b.s()); + for (src, tgt) in std::iter::zip(results, &mut check) { d.apply(tgt.as_slice_mut(), p - 1, target_b.t(), src.as_slice()); anyhow::ensure!( tgt.is_zero(), @@ -619,7 +682,7 @@ pub trait SecondaryLift: Sync + Sized { }; let mut f = save_file.create_file(dir.to_owned(), false); - for row in &results { + for row in results { row.to_bytes(&mut f).unwrap(); } drop(f); @@ -637,20 +700,24 @@ pub trait SecondaryLift: Sync + Sized { } } - Ok(homotopy + Ok(self.homotopies()[b.s()] .homotopies - .add_generators_from_rows_ooo(b.t(), results)) + .add_generators_from_rows_ooo(b.t(), results.to_vec())) + } + + /// Seed the zero homotopy at the bottom row `s = shift.s()`. The homotopies there are just zero; + /// every higher row is lifted from these. + fn seed_zero_homotopy(&self) { + let shift = self.shift(); + let h = &self.homotopies()[shift.s()]; + h.homotopies.extend_by_zero(h.composites.max_degree()); } #[tracing::instrument(skip(self))] fn compute_homotopies(&self) { let shift = self.shift(); - // When s = shift_s, the homotopies are just zero - { - let h = &self.homotopies()[shift.s()]; - h.homotopies.extend_by_zero(h.composites.max_degree()); - } + self.seed_zero_homotopy(); let min_t = self.homotopies()[shift.s()].homotopies.min_degree(); let s_range = self.homotopies().range(); @@ -659,6 +726,26 @@ pub trait SecondaryLift: Sync + Sized { sseq::coordinates::iter_s_t(&|b| self.compute_homotopy_step(b), min, max); } + /// Everything [`extend_all`](Self::extend_all) does *except* the final homotopy lift + /// ([`compute_homotopies`](Self::compute_homotopies)): initialize the homotopies, compute the + /// composites and intermediates, and seed the zero row. None of this consumes the target's + /// quasi-inverse. After calling this on each of several secondary lifts sharing a target, their + /// homotopy lifts can be batched together through + /// [`MultiLift`] (see + /// [`batch_extend_secondary`]). + fn prepare_homotopies(&self) { + self.initialize_homotopies(); + // A lift with no room to extend (its underlying reaches only its own shift filtration) has + // no secondary homotopies to compute — the zero row itself is out of range. Skip it; the + // [`SecondaryLiftable`] wrapper likewise finds no in-range step to drive. + if !self.homotopies().range().contains(&self.shift().s()) { + return; + } + self.compute_composites(); + self.compute_intermediates(); + self.seed_zero_homotopy(); + } + #[tracing::instrument(skip(self))] fn extend_all(&self) { self.initialize_homotopies(); @@ -668,6 +755,108 @@ pub trait SecondaryLift: Sync + Sized { } } +/// Outcome of [`SecondaryLift::prepare_homotopy_step`]. +pub enum SecondaryHomotopyPrep { + /// The step needed no quasi-inverse and is already finished; carries the range of + /// newly-contiguous degrees. + Done(std::ops::Range), + /// The step needs a quasi-inverse solve at `target_b`; complete it with + /// [`SecondaryLift::finish_homotopy_step`]. + NeedsLift(SecondaryPending), +} + +/// A secondary homotopy step awaiting its quasi-inverse solve; see +/// [`SecondaryLift::prepare_homotopy_step`]. +pub struct SecondaryPending { + b: Bidegree, + target_b: Bidegree, + num_gens: usize, + target_dim: usize, + /// The intermediates to lift, one per source generator. + intermediates: Vec, + /// Copy of the intermediates for the post-solve lift-validity check; `Some` only at the bottom + /// non-trivial row `b.s() == shift.s() + 1`. + check: Option>, +} + +/// A [`Liftable`] view of a secondary lift, so its homotopy solve can be driven by [`MultiLift`] +/// alongside other secondary lifts that share the same target complex. +/// +/// The bidegree `b` that [`Liftable::prepare`] receives is the *target* bidegree where the +/// quasi-inverse is taken; the secondary homotopy step it drives is at `b + shift + (0, 1)`. +struct SecondaryLiftable(Arc); + +impl Liftable for SecondaryLiftable +where + T: SecondaryLift + Send + Sync + 'static, +{ + fn prepare(&self, b: Bidegree) -> Option> { + let lift = &*self.0; + let shift = lift.shift(); + // The homotopy source bidegree whose quasi-inverse solve happens at `b`. + let source = Bidegree::s_t(b.s() + shift.s(), b.t() + shift.t() + 1); + + // Row `shift.s()` is the separately-seeded zero homotopy; outside the computed range there + // is nothing to lift. These guards keep the indexing below in bounds and stop `MultiLift` + // (which sweeps the whole target plane) from driving steps past what this lift covers. + let max = lift.max(); + if source.s() < shift.s() + 1 || source.s() >= max.s() || source.t() >= max.t(source.s()) { + return None; + } + + match lift.prepare_homotopy_step(source) { + SecondaryHomotopyPrep::Done(_) => None, + SecondaryHomotopyPrep::NeedsLift(mut pending) => { + let inputs = std::mem::take(&mut pending.intermediates); + Some(LiftRequest { + inputs, + finish: Box::new(move |results| { + // Panics on an invalid (non-lifting) input, matching the `.unwrap()` in the + // single-lift driver [`SecondaryLift::compute_homotopy_step`]. + lift.finish_homotopy_step(pending, results).unwrap(); + }), + }) + } + } + } +} + +/// Extend several secondary lifts that share one target complex together, batching the target's +/// quasi-inverse solve at each bidegree via [`MultiLift`] — the secondary analogue of +/// [`ExtAlgebra::extend_all_products`](crate::ext_algebra::ExtAlgebra::extend_all_products). +/// +/// Each lift's composites and intermediates (which do not touch the shared quasi-inverse) are +/// computed first via [`SecondaryLift::prepare_homotopies`]; only the homotopy solve is batched, so +/// with quasi-inverses recomputed on demand the target's quasi-inverse at each bidegree is solved +/// once and shared across all the lifts rather than once per lift. +/// +/// All `lifts` must lift through `target` (i.e. `lift.target()` is `target` for each). +pub fn batch_extend_secondary(target: Arc, lifts: Vec>) +where + T: SecondaryLift + Send + Sync + 'static, +{ + if lifts.is_empty() { + return; + } + // Every request is solved against `target`'s quasi-inverse, so a lift through a different target + // would be silently wrong. Enforce the shared-target contract before any preparation side + // effects, matching the `Arc::ptr_eq` checks the lift constructors already use. + for lift in &lifts { + assert!( + Arc::ptr_eq(&target, &lift.target()), + "batch_extend_secondary: every lift must share the supplied target" + ); + } + for lift in &lifts { + lift.prepare_homotopies(); + } + let liftables: Vec> = lifts + .into_iter() + .map(|l| Arc::new(SecondaryLiftable(l)) as Arc) + .collect(); + MultiLift::new(target, liftables).extend_all(); +} + #[cfg(test)] mod tests { use serde_json::json; @@ -748,4 +937,34 @@ mod tests { assert!(result.is_err()); assert!(result.unwrap_err().to_string().contains("Failed to lift")); } + + /// Driving a secondary resolution through [`batch_extend_secondary`] (i.e. via [`MultiLift`]) + /// must produce exactly the same homotopies — and hence the same $d_2$ — as the native + /// single-lift [`SecondaryLift::extend_all`]. This pins the [`SecondaryLiftable`] wrapper and its + /// bidegree mapping against the reference path. + #[test] + fn batched_matches_native_d2() { + let res = Arc::new(crate::utils::construct_standard::("S_2", None).unwrap()); + res.compute_through_stem(Bidegree::n_s(16, 6)); + + let native = SecondaryResolution::new(Arc::clone(&res)); + native.extend_all(); + + let batched = Arc::new(SecondaryResolution::new(Arc::clone(&res))); + batch_extend_secondary(Arc::clone(&res), vec![Arc::clone(&batched)]); + + // Compare the d2-defining matrix `homotopy(s + 2).hom_k(t)` at every stem where it is + // defined; equality across the plane means the two homotopies agree. + let mut compared = 0; + for b in res.iter_stem() { + if !(b.t() > 0 && res.has_computed_bidegree(b + Bidegree::n_s(-1, 2))) { + continue; + } + let native_m = native.homotopy(b.s() + 2).homotopies.hom_k(b.t()); + let batched_m = batched.homotopy(b.s() + 2).homotopies.hom_k(b.t()); + assert_eq!(native_m, batched_m, "d2 matrices differ at {b}"); + compared += 1; + } + assert!(compared > 0, "expected to compare some d2 matrices"); + } }