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
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
//! The `[n, k]` interaction layout the CPU, GPU and hardware emitters share.
//!
//! The CSR form in [`crate::graph::Graph`] is right for a CPU walking a variable-degree graph. A
//! GPU workgroup and an RTL datapath both want the opposite: a rectangle, one row per node, every
//! row the same width, so a lane's work is known before it starts. This is that rectangle.
//!
//! **On padding.** Rows are padded to the maximum degree, and the padding is described by an
//! explicit `active` mask. THRML pads with a sentinel value instead — `INVALID_BIAS = -1e10` — and
//! a sentinel is a number that is fine until someone does arithmetic on it, at which point it
//! silently dominates every real term in the sum. The mask here cannot do that, and the padding is
//! deliberately inert twice over: padded slots carry weight `0.0` *and* `active = 0`, so a kernel
//! that forgot the mask entirely still computes the right field. Defence in depth costs one `f64`
//! of zero per padded slot and removes a whole class of silent wrongness.
use crate::graph::Graph;
/// A graph as a dense `[n, k]` rectangle plus a per-slot active mask.
#[derive(Clone, Debug)]
pub struct Padded {
/// Nodes.
pub n: usize,
/// Row width: the maximum degree in the graph, so every row is `k` wide.
pub k: usize,
/// Neighbour index per slot, `n * k`. Padded slots point at themselves, which is harmless
/// because their weight is zero and their mask is clear.
pub nbr: Vec<u32>,
/// Coupling per slot, `n * k`. Padded slots are exactly `0.0`.
pub w: Vec<f64>,
/// 1 for a real coupling, 0 for padding, `n * k`.
pub active: Vec<u8>,
/// Per-node bias, `n`.
pub h: Vec<f64>,
}
impl Padded {
pub fn from_graph(g: &Graph) -> Padded {
let n = g.n;
let k = g.max_degree().max(1);
let mut nbr = vec![0u32; n * k];
let mut w = vec![0.0f64; n * k];
let mut active = vec![0u8; n * k];
for i in 0..n {
// padded slots point at their own row so an ignored mask reads a defined index
for slot in 0..k {
nbr[i * k + slot] = i as u32;
}
for (slot, kk) in (g.offset[i]..g.offset[i + 1]).enumerate() {
nbr[i * k + slot] = g.nbr[kk];
w[i * k + slot] = g.w[kk];
active[i * k + slot] = 1;
}
}
Padded { n, k, nbr, w, active, h: g.h.clone() }
}
/// Local field at node `i`: `sum_j J_ij s_j + h_i`.
pub fn field(&self, i: usize, s: &[i8]) -> f64 {
let mut f = self.h[i];
for slot in 0..self.k {
let t = i * self.k + slot;
f += self.w[t] * s[self.nbr[t] as usize] as f64;
}
f
}
/// Field computed while honouring the mask, for a kernel that branches on it.
pub fn field_masked(&self, i: usize, s: &[i8]) -> f64 {
let mut f = self.h[i];
for slot in 0..self.k {
let t = i * self.k + slot;
if self.active[t] != 0 {
f += self.w[t] * s[self.nbr[t] as usize] as f64;
}
}
f
}
/// Fraction of slots that are real. Low occupancy means a rectangle is the wrong shape for this
/// graph and the caller is paying for padding it will never use.
pub fn occupancy(&self) -> f64 {
if self.active.is_empty() {
return 1.0;
}
self.active.iter().filter(|&&a| a != 0).count() as f64 / self.active.len() as f64
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::graph::GraphBuilder;
use crate::rng::Pcg;
fn random_graph(n: usize, p: f64, seed: u64) -> Graph {
let mut rng = Pcg::new(seed, 0);
let mut b = GraphBuilder::new(n);
for i in 0..n {
for j in (i + 1)..n {
if rng.f64() < p {
b.couple(i, j, rng.f64() * 2.0 - 1.0);
}
}
b.bias(i, rng.f64() - 0.5);
}
b.build()
}
#[test]
fn the_field_matches_the_csr_form() {
// If these two disagree, the GPU and the CPU are sampling different models.
for (n, p, seed) in [(12, 0.3, 1), (40, 0.15, 2), (64, 0.5, 3), (7, 1.0, 4)] {
let g = random_graph(n, p, seed);
let d = Padded::from_graph(&g);
let mut rng = Pcg::new(seed, 99);
for _ in 0..20 {
let s: Vec<i8> = (0..n).map(|_| if rng.f64() < 0.5 { 1 } else { -1 }).collect();
for i in 0..n {
let a = g.field(i, &s);
assert!((a - d.field(i, &s)).abs() < 1e-12, "n={n} i={i}");
assert!((a - d.field_masked(i, &s)).abs() < 1e-12, "masked n={n} i={i}");
}
}
}
}
#[test]
fn padding_is_inert_even_if_the_mask_is_ignored() {
// The contrast with a -1e10 sentinel: forgetting the mask must be survivable.
let g = random_graph(30, 0.1, 7); // ragged degrees, so there is real padding
let d = Padded::from_graph(&g);
assert!(d.occupancy() < 1.0, "this test needs a graph with padding");
for t in 0..d.n * d.k {
if d.active[t] == 0 {
assert_eq!(d.w[t], 0.0, "a padded slot must weigh exactly zero");
assert!((d.nbr[t] as usize) < d.n, "a padded slot must still index in range");
}
}
// field() ignores the mask entirely and still agrees with the CSR truth
let s: Vec<i8> = (0..d.n).map(|i| if i % 3 == 0 { 1 } else { -1 }).collect();
for i in 0..d.n {
assert!((g.field(i, &s) - d.field(i, &s)).abs() < 1e-12);
}
}
#[test]
fn a_sentinel_would_have_destroyed_the_field() {
// Demonstrating why this design choice is not fussiness. Substituting THRML's sentinel into
// the padding and ignoring the mask makes the field meaningless.
let g = random_graph(30, 0.1, 7);
let mut d = Padded::from_graph(&g);
for t in 0..d.n * d.k {
if d.active[t] == 0 {
d.w[t] = -1e10;
}
}
let s: Vec<i8> = (0..d.n).map(|_| 1).collect();
let ruined = (0..d.n).any(|i| (g.field(i, &s) - d.field(i, &s)).abs() > 1.0);
assert!(ruined, "a sentinel in the padding should wreck an unmasked field, and it did");
}
#[test]
fn a_regular_graph_pads_to_nothing() {
// A lattice has uniform degree, so the rectangle is exactly the graph.
let g = crate::ising::lattice2d(8, 1.0);
let d = Padded::from_graph(&g);
assert_eq!(d.k, 4);
assert_eq!(d.occupancy(), 1.0);
}
#[test]
fn occupancy_reports_what_the_rectangle_costs() {
// One high-degree node forces every row wide; the caller deserves to know.
let mut b = GraphBuilder::new(20);
for j in 1..20 {
b.couple(0, j, 1.0); // a star
}
let d = Padded::from_graph(&b.build());
assert_eq!(d.k, 19);
assert!(d.occupancy() < 0.11, "a star should pad badly, got {}", d.occupancy());
}
}