Allow computing edge orbits in normal_form

Use this option for the graphs produced in expanding_differential, so
that in the next step, only one representative of each edge orbit can be
contracted. So far we do not use this option for the graphs produced by
contractions. (It could be used in the future to produce coboundaries.)

Cargo.lock

@@ -90,6 +90,7 @@ version = "0.1.0"
 dependencies = [
  "itertools",
  "nauty-Traces-sys",
+ "union-find",
 ]
 
 [[package]]
@@ -267,6 +268,12 @@ version = "1.0.14"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "adb9e6ca4f869e1180728b7950e35922a7fc6397f7b641499e8f3ef06e50dc83"
 
+[[package]]
+name = "union-find"
+version = "0.4.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "039142448432983c34b64739f8526f8f233a1eec7a66e61b6ab29acfa781194e"
+
 [[package]]
 name = "windows-targets"
 version = "0.52.6"

Cargo.toml

@@ -6,3 +6,4 @@ edition = "2021"
 [dependencies]
 itertools = "0.13.0"
 nauty-Traces-sys = "0.8.0"
+union-find = "0.4.3"

src/main.rs

@@ -6,9 +6,14 @@ fn main() {
     println!("g = {:?}", g);
     println!("g.to_nauty() = {:?}", g.to_nauty());
     println!("UndirectedGraph::from_nauty(&g.to_nauty()) = {:?}", UndirectedGraph::from_nauty(&g.to_nauty()));
-    let (_sign, normal_g, orbits) = g.normal_form();
+    let (_sign, normal_g, orbits,_) = g.normal_form(false);
     println!("normal_g = {:?}", normal_g);
-    println!("normal_g.expanding_differential(orbits) = {:?}", normal_g.expanding_differential(&orbits));
+    println!("normal_g.expanding_differential(&orbits) = {:?}", normal_g.expanding_differential(&orbits));
+    println!("normal_g.contractions() = {:?}", normal_g.contractions(&vec![0_usize]));
+    let h: UndirectedGraph = UndirectedGraph::new(7,vec![(0, 1), (0, 2), (0, 6), (1, 4), (1, 6), (2, 5), (2, 6), (3, 4), (3, 5), (3, 6), (4, 5)]);
+    let edge_orbit_reps_h: Vec<usize> = vec![0, 2, 3, 4, 7, 9, 10];
+    println!("normal_h.contractions(&edge_orbits_h) = {:?}", h.contractions(&edge_orbit_reps_h));
     let triangle: UndirectedGraph = UndirectedGraph::new(3, vec![(0,1), (1,2), (0,2)]);
-    println!("triangle.normal_form() = {:?}", triangle.normal_form());
+    println!("triangle.normal_form() = {:?}", triangle.normal_form(true));
+    println!("triangle.contractions() = {:?}", triangle.contractions(&vec![0_usize]));
 }

src/undirected_graph.rs

@@ -120,7 +120,10 @@ impl UndirectedGraph {
         }
     }
 
-    pub fn normal_form(&self) -> (i32, UndirectedGraph, HashMap<usize, usize>) {
+    pub fn normal_form(
+        &self,
+        compute_edge_orbits: bool,
+    ) -> (i32, UndirectedGraph, HashMap<usize, usize>, Vec<usize>) {
         let mut nauty_graph: Vec<u64> = self.to_nauty();
 
         let mut options = nauty_Traces_sys::optionblk::default();
@@ -172,27 +175,60 @@ impl UndirectedGraph {
             }
         }
 
-        let mut is_odd: bool = false;
+        let mut sign: i32 = 1;
+        let mut normal_self_edge_orbits: Vec<usize> = vec![];
         Self::AUTOMORPHSM_GROUP_GENERATORS.with_borrow(|gens| {
+            // NOTE: Use Union-Find to compute edge orbits.
+            use union_find::{QuickUnionUf, UnionByRank, UnionFind};
+            let mut uf: QuickUnionUf<UnionByRank> =
+                QuickUnionUf::<UnionByRank>::new(self.edges.len());
+
             for automorphism in gens {
                 // NOTE: The automorphisms produced by nauty act on the original graph.
                 let sigma: Vec<i32> = Self::induced_edge_permutation(&self, &self, automorphism);
+
+                if compute_edge_orbits {
+                    // Take powers of sigma and act on edges of self.
+                    let mut tau: Vec<i32> = sigma.clone();
+                    let num_edges: i32 = self.edges.len() as i32;
+                    let identity_permutation: Vec<i32> = (0..num_edges).collect();
+                    while tau != identity_permutation {
+                        for k in 0..self.edges.len() {
+                            uf.union(k, tau[k] as usize);
+                        }
+                        for e in tau.iter_mut() {
+                            *e = sigma[*e as usize];
+                        }
+                    }
+                }
+
                 if Self::permutation_sign(sigma) == -1 {
-                    is_odd = true;
+                    sign = 0;
+                    // NOTE: Giving up computation of edge orbits; we will not use them anyway.
                     break;
                 }
             }
-        });
 
-        let sign: i32 = if is_odd {
-            0
-        } else {
-            let edge_permutation: Vec<i32> =
-                Self::induced_edge_permutation(&self, &normal_self, &lab);
-            Self::permutation_sign(edge_permutation)
-        };
+            if sign != 0 {
+                let edge_permutation: Vec<i32> =
+                    Self::induced_edge_permutation(&normal_self, &self, &lab_inv);
+                if compute_edge_orbits {
+                    normal_self_edge_orbits = (0..self.edges.len()).map(|e| uf.find(e)).collect();
+                    // So far these are the edge orbits of self. Finally we apply the isomorphism to get the edge orbits of normal_self.
+                    for e in &mut normal_self_edge_orbits {
+                        *e = edge_permutation[*e] as usize;
+                    }
+                }
+                sign = Self::permutation_sign(edge_permutation);
+            }
+        });
 
-        return (sign, normal_self, normal_self_orbits);
+        return (
+            sign,
+            normal_self,
+            normal_self_orbits,
+            normal_self_edge_orbits,
+        );
     }
 
     fn expand_vertex(&self, &position: &usize) -> Vec<(usize, UndirectedGraph)> {
@@ -260,27 +296,28 @@ impl UndirectedGraph {
     pub fn expanding_differential(
         &self,
         orbits: &HashMap<usize, usize>,
-    ) -> HashMap<UndirectedGraph, i64> {
-        let mut diff: HashMap<UndirectedGraph, i64> = HashMap::new();
+    ) -> HashMap<UndirectedGraph, (i64, Vec<usize>)> {
+        let mut diff: HashMap<UndirectedGraph, (i64, Vec<usize>)> = HashMap::new();
 
         for (&orbit_rep, &orbit_size) in orbits {
             for (c, dg) in self.expand_vertex(&orbit_rep) {
                 // TODO: Avoid computing orbits when not used.
-                let (normal_dg_sign, normal_dg, _normal_dg_orbits) = dg.normal_form();
+                let (normal_dg_sign, normal_dg, _normal_dg_orbits, normal_dg_edge_orbits) =
+                    dg.normal_form(true);
                 if normal_dg_sign == 0 {
                     continue;
                 }
                 // TODO: Fix coefficient types.
                 let normal_dg_coeff: i64 =
                     (normal_dg_sign as i64) * (orbit_size as i64) * (c as i64);
-                if let Some(diff_coeff) = diff.get_mut(&normal_dg) {
+                if let Some((diff_coeff, _)) = diff.get_mut(&normal_dg) {
                     *diff_coeff += normal_dg_coeff;
                 } else {
-                    diff.insert(normal_dg, normal_dg_coeff);
+                    diff.insert(normal_dg, (normal_dg_coeff, normal_dg_edge_orbits));
                 }
             }
         }
-        diff.retain(|_dg, dg_coeff| *dg_coeff != 0);
+        diff.retain(|_dg, (dg_coeff, _)| *dg_coeff != 0);
 
         return diff;
     }
@@ -323,7 +360,7 @@ impl UndirectedGraph {
                 .all(|new_edge| unique_edges.insert(*new_edge));
             if has_no_double_edges {
                 let cg: UndirectedGraph = UndirectedGraph::new(self.num_vertices - 1, new_edges);
-                let (_normal_cg_sign, normal_cg, normal_cg_orbits) = cg.normal_form();
+                let (_normal_cg_sign, normal_cg, normal_cg_orbits, _) = cg.normal_form(false);
                 contractions.insert(normal_cg, normal_cg_orbits);
             }
         }