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);

    // Reset all forces to zero.
    for force in &mut forces[..] {
        force.reset();
    }

    // Calculate repulsive force between all pairs.
    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);
        }
    }

    // Calculate spring force between adjacent pairs.
    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 {
            // XXX: can we treat the scaled force as moved distance?
            node_positions[i].add_scaled(step / length, &forces[i]);
        }
        node_positions[i].clip_within(min_pos, max_pos);

        // add up the moved distance
        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);

    // initialize forces
    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);
}