1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
pub use vector::{Vector, P2d};
mod vector;
#[inline]
fn attractive_force<V>(p1: &V, p2: &V, k_s: f32) -> V where V: Vector<Scalar=f32>
{
let mut force = p1.sub(p2);
let length = force.length_squared().sqrt();
let strength = length / k_s;
force.scale(strength);
return force;
}
#[inline]
fn repulsive_force<V>(p1: &V, p2: &V, k_r: f32) -> V where V: Vector<Scalar=f32>
{
let mut force = p1.sub(p2);
let length_squared = force.length_squared();
if length_squared > 0.0 {
let strength = k_r / length_squared;
force.scale(strength);
}
return force;
}
fn calculate_node_forces<V>(forces: &mut[V], node_positions: &[V], node_neighbors: &[Vec<usize>],
k_r: f32, k_s: f32) where V: Vector<Scalar=f32> {
let n = forces.len();
assert!(node_positions.len() == n);
assert!(node_neighbors.len() == n);
for force in &mut forces[..] {
force.reset();
}
for i1 in 0 .. n - 1 {
for i2 in i1 + 1 .. n {
let force = repulsive_force(&node_positions[i1], &node_positions[i2], k_r);
forces[i1].add_scaled(1.0, &force);
forces[i2].add_scaled(-1.0, &force);
}
}
for i1 in 0 .. n {
for i2 in node_neighbors[i1].iter().map(|&i| i as usize) {
if i1 < i2 {
let force = attractive_force(&node_positions[i1], &node_positions[i2], k_s);
forces[i1].add_scaled(-1.0, &force);
forces[i2].add_scaled(1.0, &force);
}
}
}
}
fn update_node_positions<V>(forces: &[V], node_positions: &mut[V], step: f32, min_pos: &V, max_pos: &V) -> f32
where V: Vector<Scalar=f32> {
let n = forces.len();
assert!(node_positions.len() == n);
let mut sum_distance = 0.0;
for i in 0..n {
let before = node_positions[i].clone();
let length = forces[i].length_squared().sqrt();
if length > 0.0 {
node_positions[i].add_scaled(step / length, &forces[i]);
}
node_positions[i].clip_within(min_pos, max_pos);
sum_distance += before.sub(&node_positions[i]).length_squared().sqrt();
}
return sum_distance;
}
pub fn fruchterman_reingold<F,V>(step_fn: F, max_iter: usize, converge_eps: f32, k_r: f32, k_s: f32, min_pos: &V, max_pos: &V,
node_positions: &mut[V], node_neighbors: &[Vec<usize>])
where V: Vector<Scalar=f32>, F: Fn(usize) -> f32
{
let n = node_positions.len();
assert!(node_neighbors.len() == n);
let mut forces: Vec<V> = (0..n).map(|_| V::new()).collect();
let mut iter: usize = 0;
while iter < max_iter {
let step = step_fn(iter);
iter += 1;
calculate_node_forces(&mut forces, node_positions, node_neighbors, k_r, k_s);
let dist_moved = update_node_positions(&forces, node_positions, step, min_pos, max_pos);
if dist_moved < converge_eps {
break;
}
}
}
pub fn typical_fruchterman_reingold_2d(area: f32, node_positions: &mut[P2d], node_neighbors: &[Vec<usize>]) {
let n = node_positions.len();
assert!(node_neighbors.len() == n);
const MAX_ITER: usize = 300;
const EPS: f32 = 0.1;
const OPT_PAIR_SQR_DIST_SCALE: f32 = 0.3;
let min_pos = P2d(0.0, 0.0);
let max_pos = P2d(1.0, 1.0);
let step_fn = |iter| { 1.0 / (1.0 + (iter as f32) * 0.1) };
let k_r: f32 = OPT_PAIR_SQR_DIST_SCALE * area / (n as f32);
let k_s = k_r.sqrt();
fruchterman_reingold(step_fn, MAX_ITER, EPS, k_r, k_s, &min_pos, &max_pos, node_positions, node_neighbors);
}