From a6be55c6bc0a8e14ddde80a8490ea9efa2146f7f Mon Sep 17 00:00:00 2001 From: Antonio De Lucreziis Date: Sat, 25 May 2024 01:49:08 +0200 Subject: [PATCH] more experiments --- Cargo.lock | 15 ++ examples/fdg-example/src/main.rs | 3 + .../fdg-example/src/stress_majorization.rs | 152 +++++++++++ examples/graphs-1/Cargo.toml | 17 ++ examples/graphs-1/src/main.rs | 243 ++++++++++++++++++ src/graph.rs | 71 ++++- 6 files changed, 489 insertions(+), 12 deletions(-) create mode 100644 examples/fdg-example/src/stress_majorization.rs create mode 100644 examples/graphs-1/Cargo.toml create mode 100644 examples/graphs-1/src/main.rs diff --git a/Cargo.lock b/Cargo.lock index df86535..60e2154 100644 --- a/Cargo.lock +++ b/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" diff --git a/examples/fdg-example/src/main.rs b/examples/fdg-example/src/main.rs index ff50482..39cc4df 100644 --- a/examples/fdg-example/src/main.rs +++ b/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")] diff --git a/examples/fdg-example/src/stress_majorization.rs b/examples/fdg-example/src/stress_majorization.rs new file mode 100644 index 0000000..4c8709e --- /dev/null +++ b/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; + +use fdg::Field; +use nalgebra::Point; + +#[derive(Debug, Clone)] +pub struct StressMajorizationConfiguration { + pub dt: F, +} + +use petgraph::{ + algo::{dijkstra, Measure}, + graph::NodeIndex, + stable_graph::StableGraph, +}; +use StressMajorizationConfiguration as Config; + +impl Default for Config { + fn default() -> Self { + Self { + dt: F::from(0.035).unwrap(), + } + } +} + +/// A basic implementation +#[derive(Debug)] +pub struct StressMajorization { + pub conf: Config, + pub shortest_path_matrix: HashMap>, +} + +impl StressMajorization { + pub fn new(conf: Config) -> Self { + Self { + conf, + shortest_path_matrix: HashMap::new(), + } + } + + fn apply(&mut self, graph: &mut StableGraph<(N, Point), 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), E> = graph; + + self.shortest_path_matrix.extend( + graph + .node_indices() + .map(|idx| (idx, dijkstra(graph as &_, idx, None, |_| F::one()))), + ) + } + + fn calc_stress( + &mut self, + graph: &mut StableGraph<(N, Point), 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 Force for StressMajorization { +// fn apply(&mut self, graph: &mut ForceGraph) { +// 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, HashFn> = graph +// .node_indices() +// .map(|idx| (idx, graph.node_weight(idx).unwrap().1)) +// .collect(); + +// for idx in start_positions.keys() { +// let mut velocity: SVector = *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::>(); +// 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::>(); + +// 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); +// } +// } +// } diff --git a/examples/graphs-1/Cargo.toml b/examples/graphs-1/Cargo.toml new file mode 100644 index 0000000..80f8ef2 --- /dev/null +++ b/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" diff --git a/examples/graphs-1/src/main.rs b/examples/graphs-1/src/main.rs new file mode 100644 index 0000000..dd51f00 --- /dev/null +++ b/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::::zeros() + } + }) + .sum::>() + + { + 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::::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::::zeros() + } + }) + .sum::>(), + ) + }) + .collect::>(); + + 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), ()> { + 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), ()>) { + 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 +} diff --git a/src/graph.rs b/src/graph.rs index 876de02..73c0f4a 100644 --- a/src/graph.rs +++ b/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 { - let mut visited: HashSet = HashSet::new(); - let mut stack = vec![&node]; + pub fn dfs<'a>(&'a self, node: &'a V) -> impl Iterator + 'a { + let mut visited = HashSet::new(); + let mut stack = VecDeque::from([node]); - // O(|V| + |E|) - while let Some(node) = stack.pop() { - visited.insert(node.clone()); + 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) { - for adj in adjacencies { - if !visited.contains(adj) { - stack.push(adj); + 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; + } + + visited.insert(node.clone()); + + let distance = *distances.get(node).unwrap(); + + if let Some(adjacencies) = self.get_adjacencies(node) { + for adj in adjacencies { + if !distances.contains_key(adj) { + distances.insert(adj, distance + 1); + queue.push_back(adj); + } } } } + + result.insert(node, distances); } - visited + result } pub fn compute_ccs(&self) -> Vec> {