nassau: recompute quasi-inverses on demand, with bidegree-major batched lifting

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<u32> = 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(())
+}

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<i32> {
+        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<FpVector> = (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<i32> {
         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<FreeModuleHomomorphism<U::Module>> {
@@ -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<i32>),
+    /// 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<FpVector>,
+}
+
+impl<
+    S: FreeChainComplex + Send + Sync + 'static,
+    T: FreeChainComplex<Algebra = S::Algebra> + Send + Sync + 'static,
+    U: ChainComplex<Algebra = S::Algebra> + Send + Sync + 'static,
+> Liftable for ChainHomotopy<S, T, U>
+{
+    fn prepare(&self, b: Bidegree) -> Option<LiftRequest<'_>> {
+        // `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)`;