more experiments

Cargo.lock

@@ -1626,6 +1626,21 @@ dependencies = [
  "bitflags 2.5.0",
 ]
 
+[[package]]
+name = "graphs_1"
+version = "0.1.0"
+dependencies = [
+ "asd",
+ "fdg",
+ "macroquad",
+ "nalgebra",
+ "num-traits",
+ "petgraph",
+ "petgraph-gen",
+ "rand",
+ "rayon",
+]
+
 [[package]]
 name = "hashbrown"
 version = "0.13.2"

examples/fdg-example/src/main.rs

@@ -4,6 +4,7 @@ use asd::{
     gfa::{Entry, Orientation},
     parser,
 };
+
 use fdg::{
     fruchterman_reingold::{FruchtermanReingoldConfiguration, FruchtermanReingoldParallel},
     petgraph::Graph,
@@ -12,6 +13,8 @@ use fdg::{
     ForceGraph,
 };
 
+mod stress_majorization;
+
 use macroquad::prelude::*;
 
 #[macroquad::main("fdg demo")]

examples/fdg-example/src/stress_majorization.rs

@@ -0,0 +1,152 @@
+use std::{
+    collections::HashMap,
+    hash::{BuildHasherDefault, Hash},
+    iter::Sum,
+    ops::AddAssign,
+};
+
+// use nalgebra::{Point, SVector};
+// use petgraph::stable_graph::NodeIndex;
+
+// type HashFn = BuildHasherDefault<rustc_hash::FxHasher>;
+
+use fdg::Field;
+use nalgebra::Point;
+
+#[derive(Debug, Clone)]
+pub struct StressMajorizationConfiguration<F: Field> {
+    pub dt: F,
+}
+
+use petgraph::{
+    algo::{dijkstra, Measure},
+    graph::NodeIndex,
+    stable_graph::StableGraph,
+};
+use StressMajorizationConfiguration as Config;
+
+impl<F: Field> Default for Config<F> {
+    fn default() -> Self {
+        Self {
+            dt: F::from(0.035).unwrap(),
+        }
+    }
+}
+
+/// A basic implementation
+#[derive(Debug)]
+pub struct StressMajorization<F: Field, const D: usize> {
+    pub conf: Config<F>,
+    pub shortest_path_matrix: HashMap<NodeIndex, HashMap<NodeIndex, F>>,
+}
+
+impl<F: Field + Measure + Sum, const D: usize> StressMajorization<F, D> {
+    pub fn new(conf: Config<F>) -> Self {
+        Self {
+            conf,
+            shortest_path_matrix: HashMap::new(),
+        }
+    }
+
+    fn apply<N: Clone, E: Clone>(&mut self, graph: &mut StableGraph<(N, Point<F, D>), E>) {
+        // if self.shortest_path_matrix.is_empty() {
+        //     self.shortest_path_matrix = graph
+        //         .node_indices()
+        //         .map(|idx| {
+        //             (idx, )
+        //         })
+        //         .collect();
+        // }
+
+        // let borrowed_graph: &StableGraph<(N, Point<F, D>), E> = graph;
+
+        self.shortest_path_matrix.extend(
+            graph
+                .node_indices()
+                .map(|idx| (idx, dijkstra(graph as &_, idx, None, |_| F::one()))),
+        )
+    }
+
+    fn calc_stress<N: Clone, E: Clone>(
+        &mut self,
+        graph: &mut StableGraph<(N, Point<F, D>), E>,
+    ) -> F {
+        graph
+            .node_indices()
+            .flat_map(|v| {
+                graph.node_indices().skip(v.index() + 1).map(move |w| {
+                    let dist = nalgebra::distance(
+                        &graph.node_weight(v).unwrap().1,
+                        &graph.node_weight(w).unwrap().1,
+                    );
+
+                    if dist != F::zero() {
+                        let dij = self.shortest_path_matrix[&v][&w];
+
+                        let sp_diff = self.shortest_path_matrix[&v][&w] - dist;
+                        dij.simd_sqrt().abs() * sp_diff * sp_diff
+                    } else {
+                        F::zero()
+                    }
+                })
+            })
+            .sum()
+    }
+}
+
+// impl<const D: usize, N, E> Force<f32, D, N, E> for StressMajorization<D> {
+//     fn apply(&mut self, graph: &mut ForceGraph<f32, D, N, E>) {
+//         if self.shortest_path_matrix.is_empty() {
+//             self.shortest_path_matrix = graph
+//                 .node_indices()
+//                 .map(|idx| (idx, petgraph::algo::dijkstra(graph, idx, None, |e| 1.0)))
+//                 .collect();
+//         }
+
+//         let start_positions: HashMap<NodeIndex, Point<f32, D>, HashFn> = graph
+//             .node_indices()
+//             .map(|idx| (idx, graph.node_weight(idx).unwrap().1))
+//             .collect();
+
+//         for idx in start_positions.keys() {
+//             let mut velocity: SVector<f32, D> = *self.velocities.get(idx).unwrap_or(&SVector::<
+//                 f32,
+//                 D,
+//             >::zeros(
+//             ));
+
+//             let pos = start_positions.get(idx).unwrap();
+
+//             let attraction = graph
+//                 .neighbors_undirected(*idx)
+//                 .filter(|neighbor_idx| neighbor_idx != idx)
+//                 .map(|neighbor_idx| start_positions.get(&neighbor_idx).unwrap())
+//                 .map(|neighbor_pos| {
+//                     (neighbor_pos - pos).normalize()
+//                         * (nalgebra::distance_squared(neighbor_pos, pos) / self.conf.scale)
+//                 })
+//                 .sum::<SVector<f32, D>>();
+//             let repulsion = graph
+//                 .node_indices()
+//                 .filter(|other_idx| other_idx != idx)
+//                 .map(|other_idx| start_positions.get(&other_idx).unwrap())
+//                 .map(|other_pos| {
+//                     (other_pos - pos).normalize()
+//                         * -(self.conf.scale.simd_powi(2)
+//                             / nalgebra::distance_squared(other_pos, pos))
+//                 })
+//                 .sum::<SVector<f32, D>>();
+
+//             velocity.add_assign((attraction + repulsion) * self.conf.dt);
+//             velocity.scale_mut(self.conf.cooloff_factor);
+
+//             self.velocities.insert(*idx, velocity);
+
+//             graph
+//                 .node_weight_mut(*idx)
+//                 .unwrap()
+//                 .1
+//                 .add_assign(velocity * self.conf.dt);
+//         }
+//     }
+// }

examples/graphs-1/Cargo.toml

@@ -0,0 +1,17 @@
+[package]
+name = "graphs_1"
+version = "0.1.0"
+edition = "2021"
+
+[dependencies]
+nalgebra = { version = "0.32.3", features = ["rand"] }
+petgraph = { version = "0.6.4", features = [
+    "stable_graph",
+], default-features = false }
+num-traits = "0.2.17"
+rand = "0.8.5"
+macroquad = "0.4.4"
+petgraph-gen = "0.1.3"
+fdg = { git = "https://github.com/grantshandy/fdg" }
+asd = { path = "../../" }
+rayon = "1.10.0"

examples/graphs-1/src/main.rs

@@ -0,0 +1,243 @@
+use std::{collections::HashMap, env, ops::AddAssign, time::Instant};
+
+use asd::{gfa::Entry, parser};
+use macroquad::{prelude::*, rand};
+use nalgebra::{Point2, SVector};
+use petgraph::{algo::dijkstra, stable_graph::StableGraph};
+
+use rayon::prelude::*;
+
+#[macroquad::main("graphs_1")]
+async fn main() {
+    println!("Hello, world!");
+
+    let mut graph = load_graph();
+
+    let mut desired_distance_matrix = HashMap::new();
+
+    graph.node_indices().for_each(|idx| {
+        desired_distance_matrix.insert(idx, dijkstra(&graph, idx, None, |_| 1.0 as f32));
+    });
+
+    // println!("{:?}", desired_distance_matrix);
+
+    loop {
+        // Update
+
+        for _ in 1..2 {
+            let forces = graph
+                .node_indices()
+                .par_bridge()
+                .map(|idx| {
+                    (
+                        idx,
+                        desired_distance_matrix
+                            .get(&idx)
+                            .unwrap()
+                            .par_iter()
+                            .map(|(&other_idx, &distance)| {
+                                let pos = graph.node_weight(idx).unwrap().1;
+                                let other_pos = graph.node_weight(other_idx).unwrap().1;
+                                let delta = other_pos - pos;
+                                let dist = delta.norm_squared();
+                                let correction = dist - (distance * distance);
+
+                                // println!("correction: {:?}", correction);
+
+                                if distance > 0.0 && dist > 1e-6 {
+                                    0.001 * delta.normalize() * correction.clamp(-10.0, 10.0)
+                                } else {
+                                    SVector::<f32, 2>::zeros()
+                                }
+                            })
+                            .sum::<SVector<f32, 2>>()
+                            + {
+                                let pos = graph.node_weight(idx).unwrap().1;
+
+                                let delta = pos - Point2::new(0.0, 0.0);
+
+                                let dist = delta.norm_squared();
+
+                                if dist > 1e-6 {
+                                    -0.01 * delta.normalize() / dist.max(1.0)
+                                } else {
+                                    SVector::<f32, 2>::zeros()
+                                }
+                            }
+                            + graph
+                                .node_indices()
+                                .par_bridge()
+                                .filter(|&other_idx| other_idx != idx)
+                                .map(|other_idx| {
+                                    let pos = graph.node_weight(idx).unwrap().1;
+                                    let other_pos = graph.node_weight(other_idx).unwrap().1;
+                                    let delta = other_pos - pos;
+                                    let dist = delta.norm();
+                                    if dist > 1e-3 {
+                                        -0.001 * delta.normalize() / (dist * dist)
+                                    } else {
+                                        SVector::<f32, 2>::zeros()
+                                    }
+                                })
+                                .sum::<SVector<f32, 2>>(),
+                    )
+                })
+                .collect::<HashMap<_, _>>();
+
+            forces.iter().for_each(|(idx, force)| {
+                let (_, pos) = graph.node_weight_mut(*idx).unwrap();
+                pos.add_assign(force);
+            });
+        }
+
+        // Render
+
+        let now = Instant::now();
+
+        clear_background(WHITE);
+
+        draw_graph(&graph);
+
+        let elapsed = now.elapsed();
+
+        println!("frame: {:?}", elapsed);
+
+        next_frame().await
+    }
+}
+
+fn load_graph() -> StableGraph<(String, Point2<f32>), ()> {
+    let mut graph = StableGraph::new();
+
+    let file = std::fs::File::open(env::args().nth(1).expect("missing gfa file argument")).unwrap();
+    let entries = parser::parse_source(file).unwrap();
+
+    let mut index_map = HashMap::new();
+
+    // let node_count = entries
+    //     .iter()
+    //     .filter_map(|entry| match entry {
+    //         Entry::Segment { id, .. } => Some(id),
+    //         _ => None,
+    //     })
+    //     .count();
+
+    // let radius = (node_count as f32).sqrt();
+
+    let mut i = -10.0;
+
+    for entry in entries
+        .iter()
+        .filter(|entry| matches!(entry, Entry::Link { .. }))
+        .take(3000)
+    {
+        // println!("{:?}", entry);
+
+        if let Entry::Link {
+            from,
+            from_orient,
+            to,
+            to_orient,
+        } = entry
+        {
+            // add first node if not present
+            let a = index_map
+                .entry(from.clone())
+                .or_insert_with(|| {
+                    i += 1.0;
+
+                    graph.add_node((
+                        format!("{}{}", from, from_orient),
+                        // Point2::new(rand::gen_range(0.0, radius), rand::gen_range(0.0, radius)),
+                        Point2::new(i, 50.0 + rand::gen_range(0.0, 100.0)),
+                    ))
+                })
+                .to_owned();
+
+            // add second node if not present
+            let b = index_map
+                .entry(to.clone())
+                .or_insert_with(|| {
+                    i += 1.0;
+
+                    graph.add_node((
+                        format!("{}{}", from, to_orient),
+                        // Point2::new(rand::gen_range(0.0, radius), rand::gen_range(0.0, radius)),
+                        Point2::new(i, 50.0 + rand::gen_range(0.0, 100.0)),
+                    ))
+                })
+                .to_owned();
+
+            graph.add_edge(a, b, ());
+            graph.add_edge(b, a, ());
+        }
+    }
+
+    graph
+}
+
+fn draw_graph(graph: &StableGraph<(String, Point2<f32>), ()>) {
+    let (width, height) = (screen_width(), screen_height());
+
+    let (min_x, max_x) = graph
+        .node_weights()
+        .map(|(_, pos)| pos.x)
+        .fold((f32::INFINITY, f32::NEG_INFINITY), |(min, max), x| {
+            (min.min(x), max.max(x))
+        });
+
+    let (min_y, max_y) = graph
+        .node_weights()
+        .map(|(_, pos)| pos.y)
+        .fold((f32::INFINITY, f32::NEG_INFINITY), |(min, max), y| {
+            (min.min(y), max.max(y))
+        });
+
+    let source_range: f32 = (max_x - min_x).max(max_y - min_y);
+
+    for idx in graph.edge_indices() {
+        let ((_, source), (_, target)) = graph
+            .edge_endpoints(idx)
+            .map(|(a, b)| (graph.node_weight(a).unwrap(), graph.node_weight(b).unwrap()))
+            .unwrap();
+
+        draw_line(
+            remap(source.x, min_x, min_x + source_range, 10.0, width - 10.0),
+            remap(source.y, min_y, min_y + source_range, 10.0, height - 10.0),
+            remap(target.x, min_x, min_x + source_range, 10.0, width - 10.0),
+            remap(target.y, min_y, min_y + source_range, 10.0, height - 10.0),
+            1.0,
+            BLACK,
+        );
+    }
+
+    for (_label, pos) in graph.node_weights() {
+        let x = remap(pos.x, min_x, min_x + source_range, 10.0, width - 10.0);
+        let y = remap(pos.y, min_y, min_y + source_range, 10.0, height - 10.0);
+
+        draw_circle(x, y, 2.0, RED);
+        // draw_text(label.as_str(), x - 30.0, y - 30.0, 10.0, BLACK);
+    }
+
+    draw_line(
+        remap(0.0, min_x, min_x + source_range, 10.0, width - 10.0),
+        remap(0.0, min_y, min_y + source_range, 10.0, height - 10.0),
+        remap(100.0, min_x, min_x + source_range, 10.0, width - 10.0),
+        remap(0.0, min_y, min_y + source_range, 10.0, height - 10.0),
+        2.0,
+        BLUE,
+    );
+
+    draw_line(
+        remap(0.0, min_x, min_x + source_range, 10.0, width - 10.0),
+        remap(0.0, min_y, min_y + source_range, 10.0, height - 10.0),
+        remap(0.0, min_x, min_x + source_range, 10.0, width - 10.0),
+        remap(100.0, min_y, min_y + source_range, 10.0, height - 10.0),
+        2.0,
+        BLUE,
+    );
+}
+
+fn remap(value: f32, from_min: f32, from_max: f32, to_min: f32, to_max: f32) -> f32 {
+    (value - from_min) / (from_max - from_min) * (to_max - to_min) + to_min
+}

src/graph.rs

@@ -1,6 +1,6 @@
 use std::{
     cell::RefCell,
-    collections::{HashMap, HashSet},
+    collections::{HashMap, HashSet, VecDeque},
     fmt::Debug,
     hash::Hash,
     rc::Rc,
@@ -69,6 +69,18 @@ where
         opposite
     }
 
+    pub fn undirected(&self) -> AdjacencyGraph<&V> {
+        let mut undirected = AdjacencyGraph::new();
+
+        // O(|E|)
+        for (from, to) in self.edges() {
+            undirected.add_edge(from, to);
+            undirected.add_edge(to, from);
+        }
+
+        undirected
+    }
+
     pub fn has_edge(&self, from: &V, to: &V) -> bool {
         // O(1)
         if let Some(adjacencies) = self.get_adjacencies(from) {
@@ -79,24 +91,59 @@ where
         }
     }
 
-    pub fn dfs(&self, node: V) -> HashSet<V> {
+    pub fn dfs<'a>(&'a self, node: &'a V) -> impl Iterator<Item = V> + 'a {
-        let mut visited: HashSet<V> = HashSet::new();
+        let mut visited = HashSet::new();
-        let mut stack = vec![&node];
+        let mut stack = VecDeque::from([node]);
+
+        std::iter::from_fn(move || {
+            while let Some(node) = stack.pop_back() {
+                if !visited.insert(node.clone()) {
+                    continue;
+                }
+
+                if let Some(adjacencies) = self.get_adjacencies(node) {
+                    stack.extend(adjacencies);
+                }
+
+                return Some(node.clone());
+            }
+            None
+        })
+    }
+
+    pub fn shortest_path_matrix(&self) -> HashMap<&V, HashMap<&V, usize>> {
+        let mut result = HashMap::new();
+
+        for node in self.nodes.iter() {
+            let mut distances = HashMap::new();
+            let mut visited = HashSet::new();
+            let mut queue = VecDeque::from([node]);
+
+            distances.insert(node, 0);
+
+            while let Some(node) = queue.pop_front() {
+                if visited.contains(node) {
+                    continue;
+                }
 
-        // O(|V| + |E|)
-        while let Some(node) = stack.pop() {
                 visited.insert(node.clone());
 
-            if let Some(adjacencies) = self.get_adjacencies(&node) {
+                let distance = *distances.get(node).unwrap();
+
+                if let Some(adjacencies) = self.get_adjacencies(node) {
                     for adj in adjacencies {
-                    if !visited.contains(adj) {
+                        if !distances.contains_key(adj) {
-                        stack.push(adj);
+                            distances.insert(adj, distance + 1);
+                            queue.push_back(adj);
+                        }
                     }
                 }
             }
+
+            result.insert(node, distances);
         }
 
-        visited
+        result
     }
 
     pub fn compute_ccs(&self) -> Vec<Vec<V>> {