explaining affinity propagation
How exemplars emerge
Affinity propagation (Frey and Dueck, Science 2007, pdf) is a clustering algorithm based on message passing that magically finds a "natural" number of clusters in a dataset.
Every point exchanges two kinds of messages (responsibilities, availabilities) with every other point until a few of them stand out as exemplars (cluster centers).
This is an interactive explanation of how exemplars emerge, inspired by the original paper's Fig. 1A.
Press Play, step with the arrows (or the arrow keys), or drag the slider.
evidence r(k,k) + a(k,k) that k is an exemplar (it is one when positive)i → k: i's belief that k is its exemplar
parameters preference median similarity · damping λ = 0.9 · convits 10
The code
From the algorithm to an interactive explanation. The source of this page itself is under Show Source at the bottom.
The algorithm. naffprop's own affinity propagation, the code that runs when you
call it from Python. Each iteration updates the responsibilities r, then the
availabilities a; a point is an exemplar while r(k,k) + a(k,k) > 0. The onIteration
hook, compiled only when given, is how this page records the messages.
src/naffprop/ap.nim · 275 lines
## Affinity propagation (Frey and Dueck, Science 2007) on a dense similarity
## matrix.
##
## The algorithm is the one of Frey and Dueck's reference MATLAB code, which
## both R's apcluster and scikit-learn follow: the same message updates,
## damping, convergence test and final refinement of the exemplars. The
## preference range is the one of Frey and Dueck's preferenceRange.m.
##
## Pure Nim (no nimpy): `core.nim` exposes it to Python.
import std/[fenv, math]
import matrix, rng
export matrix
template maxFloat(T): untyped = maximumPositiveValue(T)
type
ApResult* = object
exemplars*: seq[int] ## indices of the exemplars, increasing
labels*: seq[int] ## for each point, the position of its exemplar in
## `exemplars`; all -1 if there are no exemplars
iterations*: int
converged*: bool
dpsim*: float ## sum of the similarities of points to their exemplar
expref*: float ## sum of the preferences of the exemplars
netsim*: float ## dpsim + expref, the objective AP maximizes
type Convergence* = object
## Which points were exemplars in the last `convits` iterations (a ring
## buffer) and how many times: converged when no point changed its status
## in that window.
n, convits: int
hist: seq[uint8]
count: seq[int]
k*: int ## exemplars in the last iteration
unconverged*: bool
proc initConvergence*(n, convits: int): Convergence =
Convergence(n: n, convits: convits, hist: newSeq[uint8](n * convits),
count: newSeq[int](n))
proc update*(c: var Convergence, it: int, isEx: openArray[bool], maxits: int): bool =
## Records the exemplars of iteration `it` (from 0). True when it is time to
## stop: converged with at least one exemplar, or `maxits` iterations.
c.unconverged = false
c.k = 0
let slot = it mod c.convits
for i in 0 ..< c.n:
let ex = uint8(isEx[i])
c.count[i] += ex.int - c.hist[i * c.convits + slot].int
c.hist[i * c.convits + slot] = ex
if c.count[i] > 0 and c.count[i] < c.convits: c.unconverged = true
c.k += ex.int
if it >= c.convits - 1 or it >= maxits - 1:
return (not c.unconverged and c.k > 0) or it >= maxits - 1
proc toResult*(c: Convergence, it: int): ApResult =
## The result after iteration `it`, before the exemplars are refined: all
## labels -1 and NaN similarities if there are no exemplars.
result.iterations = it + 1
result.converged = not c.unconverged and c.k > 0
result.labels = newSeq[int](c.n)
if c.k == 0:
for l in result.labels.mitems: l = -1
result.dpsim = NaN
result.expref = NaN
result.netsim = NaN
proc addNoise*[T](s: MatrixView[T], seed: int64) =
## Adds a tiny gaussian noise to every similarity, as R and scikit-learn do,
## so that ties (e.g. duplicated points) don't make the messages oscillate.
## Its scale is the precision of T, so that float32 similarities get it too.
var r = initRng(seed)
for k in 0 ..< s.n * s.n:
let x = s.data[k]
s.data[k] = x + T((epsilon(T) * x + minimumPositiveValue(T) * 100) * r.gauss())
proc setPreferences*[T](s: MatrixView[T], p: openArray[float]) =
## Puts the preferences on the diagonal: one for every point, or the same
## for all if `p` has a single value.
for i in 0 ..< s.n:
s[i, i] = T(if p.len == 1: p[0] else: p[i])
proc prepare*[T](s: MatrixView[T], p: openArray[float], noise: bool, seed: int64) =
## Readies the similarities for `affinityPropagation`, in place: noise,
## preferences on the diagonal, -Inf and NaN to -maxFloat (as R does).
if noise: s.addNoise(seed)
s.setPreferences(p)
for k in 0 ..< s.n * s.n:
let x = s.data[k]
if x.isNaN or x < -maxFloat(T): s.data[k] = -maxFloat(T)
{.push checks: off.} # hot loops: indices are in range by construction
proc finish*[T](s: MatrixView[T], conv: Convergence, it: int,
isEx: var seq[bool]): ApResult =
## The result once the messages have converged, with `isEx` the exemplars
## of the last iteration: each point to its most similar exemplar, then the
## exemplars refined and the points assigned again, as apcluster.m does.
result = conv.toResult(it)
if conv.k == 0: return
let n = s.n
let k = conv.k
var ex = newSeqOfCap[int](k)
for i in 0 ..< n:
if isEx[i]: ex.add i
proc assign(s: MatrixView[T], ex: seq[int], c: var seq[int]) =
## Each point to its most similar exemplar (the first on ties), each
## exemplar to itself.
for i in 0 ..< s.n:
var best = T(-Inf)
c[i] = 0
for e, j in ex:
if s[i, j] > best:
best = s[i, j]
c[i] = e
for e, j in ex: c[j] = e
var c = newSeq[int](n)
assign(s, ex, c)
# Refinement: in each cluster, the exemplar becomes the member with the
# largest sum of similarities from the other members (and its preference).
var members = newSeq[seq[int]](k)
for i in 0 ..< n: members[c[i]].add i
for e in 0 ..< k:
var best = -Inf
for j in members[e]:
var t = 0.0
for i in members[e]: t += float(s[i, j])
if t > best:
best = t
ex[e] = j
assign(s, ex, c)
# Exemplars in increasing order, labels pointing into them.
for j in 0 ..< n: isEx[j] = false
for j in ex: isEx[j] = true
var pos = newSeq[int](n)
for j in 0 ..< n:
if isEx[j]:
pos[j] = result.exemplars.len
result.exemplars.add j
for i in 0 ..< n:
let j = ex[c[i]]
result.labels[i] = pos[j]
if i == j: result.expref += float(s[i, i])
else: result.dpsim += float(s[i, j])
result.netsim = result.dpsim + result.expref
type OnIteration*[T] = proc (it: int, r, a: Matrix[T], isEx: seq[bool]) {.closure.}
## Called after each iteration `it` (from 0) with the responsibilities, the
## availabilities and which points are exemplars (r(k, k) + a(k, k) > 0):
## to record the messages, e.g. for the interactive explanation in docs/.
proc affinityPropagationImpl[T; record: static bool](s: MatrixView[T],
maxits, convits: int, damping: float, onIteration: OnIteration[T]): ApResult =
## The algorithm of `affinityPropagation`. `record` is known at compile
## time: without it the call to `onIteration` is not even compiled, so the
## plain run is the same code as if the hook did not exist.
let n = s.n
let lam = T(damping)
var
r = zeros[T](n) # responsibilities
a = zeros[T](n) # availabilities
colsum = newSeq[float](n)
conv = initConvergence(n, convits)
isEx = newSeq[bool](n)
it = 0
while true:
# Responsibilities, by row: r(i, j) = s(i, j) - max over j' != j of
# (a(i, j') + s(i, j')). On the way, the column sums of the positive
# responsibilities (and the diagonal) that the availabilities need.
for j in 0 ..< n: colsum[j] = 0
for i in 0 ..< n:
var max1, max2 = T(-Inf)
var jmax = 0
for j in 0 ..< n:
let v = a[i, j] + s[i, j]
if v > max1:
max2 = max1
max1 = v
jmax = j
elif v > max2:
max2 = v
for j in 0 ..< n:
let v = (1 - lam) * (s[i, j] - (if j == jmax: max2 else: max1)) +
lam * r[i, j]
r[i, j] = min(v, maxFloat(T))
if r[i, j] > 0 or i == j: colsum[j] += r[i, j]
# Availabilities: a(i, j) = min(0, r(j, j) + sum of the positive r(i', j)
# for i' not in {i, j}), and a(j, j) = sum of the positive r(i', j), i' != j.
for i in 0 ..< n:
for j in 0 ..< n:
let v = r[i, j]
var x = colsum[j] - (if v > 0 or i == j: float(v) else: 0.0)
if x > 0 and i != j: x = 0
a[i, j] = (1 - lam) * T(x) + lam * a[i, j]
# Exemplars: a(i, i) + r(i, i) > 0. Converged when no point changed its
# status in the last convits iterations.
for i in 0 ..< n: isEx[i] = a[i, i] + r[i, i] > 0
when record: onIteration(it, r, a, isEx)
if conv.update(it, isEx, maxits): break
inc it
result = finish(s, conv, it, isEx)
proc affinityPropagation*[T](s: MatrixView[T], maxits = 1000, convits = 100,
damping = 0.9): ApResult {.inline.} =
## Affinity propagation on the similarities `s`, whose diagonal holds the
## preferences. Stops when the exemplars have not changed for `convits`
## iterations, or after `maxits` iterations. Similarities must be finite:
## `prepare` turns -Inf (a pair never to link) into -maxFloat.
## The messages are stored in T (float32 halves the memory); the column
## sums are accumulated in float64.
affinityPropagationImpl[T, false](s, maxits, convits, damping, nil)
proc affinityPropagation*[T](s: MatrixView[T], maxits, convits: int,
damping: float, onIteration: OnIteration[T]): ApResult =
## The same, calling `onIteration` with the messages after every iteration.
affinityPropagationImpl[T, true](s, maxits, convits, damping, onIteration)
proc preferenceRange*[T](s: MatrixView[T], exact = false): (float, float) =
## The preferences between which AP finds from 1 or 2 clusters (the lower
## bound) to n clusters (the upper bound), ignoring the diagonal of `s`.
## The lower bound is exact with `exact = true` (O(n^3)), otherwise it is a
## cheaper, smaller bound (O(n^2)), as in Frey and Dueck's preferenceRange.m.
## -Inf similarities are skipped in the sums.
let n = s.n
template sim(i, j: int): float = (if i == j: 0.0 else: float(s[i, j]))
template addFinite(acc: var float, v: float) =
if v > -Inf:
acc = (if acc == -Inf: v else: acc + v)
var pmax = -Inf
var dpsim1 = -Inf # best net similarity with one exemplar
for j in 0 ..< n:
var t = -Inf
for i in 0 ..< n:
t.addFinite sim(i, j)
if i != j and s[i, j] > pmax: pmax = float(s[i, j])
if t > dpsim1: dpsim1 = t
var pmin: float
if dpsim1 == -Inf:
pmin = NaN
elif exact:
var dpsim2 = -Inf # best net similarity with two exemplars
for j1 in 0 ..< n - 1:
for j2 in j1 + 1 ..< n:
var t = -Inf
for i in 0 ..< n:
t.addFinite max(sim(i, j1), sim(i, j2))
if t > dpsim2: dpsim2 = t
pmin = dpsim1 - dpsim2
else:
# Upper bound on dpsim2: every point to its most similar other point,
# except the two for which that is least similar (the two exemplars).
var total = -Inf
var m1, m2 = Inf
for i in 0 ..< n:
var m = -Inf
for j in 0 ..< n:
if j != i and s[i, j] > m: m = float(s[i, j])
if m > -Inf:
total.addFinite m
if m < m1:
m2 = m1
m1 = m
elif m < m2:
m2 = m
pmin = if m2 == Inf: -Inf else: dpsim1 - total + m1 + m2
(pmin, pmax)
{.pop.}
The data. The types of the explanation (points, parameters, the messages at every
iteration) and the run behind the figure: the paper's own 25 points (toy25.nim),
similarities -‖xᵢ − xₖ‖², and AP with the hook recording every iteration. It runs beforehand, in C; the page reads the
result as JSON.
The points are the authors' ToyProblemData.txt, from jincheng9/AffinityPropagation.
data.nim · 130 lines
## Types and data of the interactive explanation of affinity propagation (the
## paper's Fig. 1: points in the plane, the messages at every iteration).
##
## All the data is computed beforehand, in C, by naffprop's own `ap.nim` (its
## `onIteration` hook records the messages), and written as JSON; the page only
## reads it. The types compile to JS too, so the page parses the JSON with
## std/json's `to`.
import std/[algorithm, math]
import matrix, ap
export matrix, ap
type
Point* = object
x*, y*: float
Ap2DInput* = object
## Points in the plane and their similarities.
points*: seq[Point]
similarityMatrix*: Matrix[float] ## negative squared distances, 0 on the
## diagonal: the preference goes there
## only when AP runs
ApParameters* = object
damping*: float ## lambda: each message is damping * old + (1 - damping) * new
quantile*: float ## preference = this quantile of the similarities
## between different points (0.5: the median, the paper's)
maxits*: int ## at most this many iterations
convits*: int ## stop when the exemplars stay the same this many iterations
ApIterations* = object
## The messages after each iteration (from the first: before it, all are
## 0), and the decisions they imply at that time.
r*, a*: seq[Matrix[float]] ## responsibilities r(i, k), availabilities a(i, k)
exemplars*: seq[seq[int]] ## the points k with r(k, k) + a(k, k) > 0
choice*: seq[seq[int]] ## for each point i, the k maximizing
## a(i, k) + r(i, k): i itself if it is an
## exemplar, else the exemplar it picks
ApRun* = object
## Everything a page needs about one run.
input*: Ap2DInput
params*: ApParameters
preference*: float ## the value the quantile gives, for every point
iterations*: ApIterations
result*: ApResult ## naffprop's result: the exemplars after the last
## iteration, refined (so they can differ from the
## last `iterations.exemplars`, which are not)
func paperParameters*(): ApParameters =
## The paper's Fig. 1: damping 0.5, the median similarity as preference;
## maxits and convits as naffprop's (and R's) defaults.
ApParameters(damping: 0.5, quantile: 0.5, maxits: 1000, convits: 100)
func fig1Parameters*(): ApParameters =
## The animation's: more damping than the paper's, so that the exemplars
## emerge over more iterations, and convits 10, so that it stops soon
## after (with 0.9, new exemplars appear 6 or 7 iterations apart: fewer
## stops before the third one).
ApParameters(damping: 0.9, quantile: 0.5, maxits: 1000, convits: 10)
func negDistMatrix*(points: seq[Point]): Matrix[float] =
## Negative squared Euclidean distances, the similarity of the paper's Fig. 1.
let n = points.len
result = zeros[float](n)
for i in 0 ..< n:
for k in 0 ..< n:
result[i, k] = -((points[i].x - points[k].x) ^ 2 + (points[i].y - points[k].y) ^ 2)
func initAp2DInput*(points: seq[Point]): Ap2DInput =
Ap2DInput(points: points, similarityMatrix: negDistMatrix(points))
func quantile*(xs: seq[float], q: float): float =
## As numpy's default (linear interpolation), so as naffprop's Python side.
let v = sorted(xs)
let h = q * float(v.len - 1)
let lo = int(floor(h))
let hi = min(lo + 1, v.len - 1)
v[lo] + (h - float(lo)) * (v[hi] - v[lo])
func preferenceOf*(input: Ap2DInput, params: ApParameters): float =
## The quantile `params.quantile` of the similarities between different points.
let s = input.similarityMatrix
var xs: seq[float]
for i in 0 ..< s.n:
for k in 0 ..< s.n:
if i != k: xs.add s[i, k]
quantile(xs, params.quantile)
proc runAp*(input: Ap2DInput, params: ApParameters): ApRun =
## Runs naffprop's AP (no noise: the data is deterministic) and records the
## messages and decisions after every iteration.
result = ApRun(input: input, params: params, preference: preferenceOf(input, params))
var s = input.similarityMatrix
let v = s.view
v.prepare([result.preference], noise = false, seed = 0)
var its: ApIterations
proc record(it: int, r, a: Matrix[float], isEx: seq[bool]) =
its.r.add r
its.a.add a
var ex: seq[int]
var choice = newSeq[int](r.n)
for i in 0 ..< r.n:
if isEx[i]: ex.add i
var best = -Inf
for k in 0 ..< r.n:
if a[i, k] + r[i, k] > best:
best = a[i, k] + r[i, k]
choice[i] = k
its.exemplars.add ex
its.choice.add choice
result.result = v.affinityPropagation(params.maxits, params.convits,
params.damping, record)
result.iterations = its
func roundSignificant*(x: float, digits: int): float =
## x to `digits` significant digits (the messages don't need more on a page,
## and the JSON is much smaller).
if x == 0 or x.isNaN or x.classify in {fcInf, fcNegInf}: return x
let e = floor(log10(abs(x)))
let f = pow(10.0, float(digits - 1) - e)
round(x * f) / f
proc roundMessages*(run: var ApRun, digits = 4) =
for m in run.iterations.r.mitems:
for x in m.data.mitems: x = roundSignificant(x, digits)
for m in run.iterations.a.mitems:
for x in m.data.mitems: x = roundSignificant(x, digits)
toy25.nim · 27 lines
## The paper's 25 points (Frey and Dueck 2007, Fig. 1): the toy data set the
## explanation runs on. These are the original ToyProblemData.txt values from
## the authors' AP web page (psi.toronto.edu, now offline), as mirrored in
## https://github.com/jincheng9/AffinityPropagation (commit ad82f68, 2014).
## An affine fit onto the paper's Fig. 1A lands every point within 4 px.
## With negative squared distances and the median preference, AP picks points
## 2, 6 and 19 (from 0) as exemplars, as in the paper.
import data
const toy25Points* = [
(-2.341500, 3.696800), (-1.109200, 3.111700), (-1.566900, 1.835100), (-2.658500, 0.664900),
(-4.031700, 2.845700), (-3.081000, 2.101100), (2.588000, 1.781900), (3.292300, 3.058500),
(4.031700, 1.622300), (3.081000, -0.611700), (0.264100, 0.398900), (1.320400, 2.207400),
(0.193700, 3.643600), (1.954200, -0.505300), (1.637300, 1.409600), (-0.123200, -1.516000),
(-1.355600, -3.058500), (0.017600, -4.016000), (1.003500, -3.590400), (0.017600, -2.420200),
(-1.531700, -0.930900), (-1.144400, 0.505300), (0.616200, -1.516000), (1.707700, -2.207400),
(2.095100, 3.430900)]
const paperExemplars* = [2, 6, 19]
func toy25*(): Ap2DInput =
## The points, with negative squared distances as similarities.
var points: seq[Point]
for (x, y) in toy25Points: points.add Point(x: x, y: y)
initAp2DInput(points)
The animation. A Karax app compiled to JavaScript: from the recorded messages it draws the points, the arrows and the color bar of one iteration, and the controls move through them. The figure around it and its legend are plain HTML, written at build time.
animation.nim · 263 lines
## The paper's Fig. 1A as an animation: the points of an `ApRun`,
## colored by the evidence that each is an exemplar, with an arrow from i to k
## as dark as i's belief that k is its exemplar, and under them a color bar of
## the evidence with every point's current value; play, step or scrub through
## the iterations. A Karax app (JS only) mounted on a div of the page:
## `figureHtml(id, params)` writes the div and its caption at build time,
## `mountAnimation(id, run)` starts the app in the browser.
##
## The evidence colors are fixed (they read on light and dark backgrounds);
## the rest come from CSS custom properties (--ap-ink, --ap-frame, --muted),
## so the page decides them, for both themes.
import std/[math, strutils]
import data
when defined(js):
include karax/prelude
import karax/[kdom, vstyles]
{.warning[CStringConv]: off.} # attribute values are cstrings in JS
func f(x: float): string = formatFloat(x, ffDecimal, 2)
type Frame = object
## The square of data coordinates shown, and its size in SVG units.
x0, y0, side, px: float
func frameOf(points: seq[Point], px: float): Frame =
var lo = Inf
var hi = -Inf
for p in points:
lo = min(lo, min(p.x, p.y))
hi = max(hi, max(p.x, p.y))
let pad = 0.08 * (hi - lo)
Frame(x0: lo - pad, y0: lo - pad, side: hi - lo + 2 * pad, px: px)
func sx(fr: Frame, x: float): float = (x - fr.x0) / fr.side * fr.px
func sy(fr: Frame, y: float): float = fr.px - (y - fr.y0) / fr.side * fr.px
# Evidence colors. Negative evidence follows viridis, from blue (the lowest
# evidence of the run) through green to yellow (just below 0): a perceptually
# uniform scale with no red in it. Positive evidence, an exemplar, is red, from
# a light red at 0 to a deep red at the highest evidence of the run.
const viridis = [(59, 82, 139), (44, 114, 142), (33, 145, 140), (39, 173, 129),
(94, 201, 98), (170, 220, 50), (253, 231, 37)] # 0.25 .. 1
const redLow = (255, 107, 94)
const redHigh = (178, 24, 24)
func mix(a, b: (int, int, int), t: float): string =
func c(x, y: int): int = int(round(float(x) + (float(y) - float(x)) * t))
"rgb(" & $c(a[0], b[0]) & "," & $c(a[1], b[1]) & "," & $c(a[2], b[2]) & ")"
func evidenceColor*(e, lo, hi: float): string =
## The color of evidence `e` on a run whose evidence spans lo < 0 < hi.
if e > 0:
mix(redLow, redHigh, clamp(e / hi, 0.0, 1.0))
else:
let x = clamp(1 - e / lo, 0.0, 1.0) * float(viridis.high)
let i = min(int(x), viridis.high - 1)
mix(viridis[i], viridis[i + 1], x - float(i))
func evidence(r, a: Matrix[float], k: int): float = r[k, k] + a[k, k]
func beliefs*(r, a: Matrix[float], i: int, temperature: float): seq[float] =
## How strongly point i believes each k is its exemplar, summing to 1: the
## messages are log-probabilities (up to a constant), so a softmax of
## a(i, k) + r(i, k) on the scale `temperature`. Before any message, every
## candidate gets 1/n; at convergence, the chosen exemplar gets almost 1.
let n = r.n
result = newSeq[float](n)
var best = -Inf
for k in 0 ..< n: best = max(best, a[i, k] + r[i, k])
var t = 0.0
for k in 0 ..< n:
result[k] = exp((a[i, k] + r[i, k] - best) / temperature)
t += result[k]
for k in 0 ..< n: result[k] /= t
# The figure's frame and caption, in HTML at build time: the legend, then
# the parameters of the run.
const evidenceDot* = """<svg width="16" height="16" viewBox="0 0 16 16" aria-hidden="true"><defs><linearGradient id="ap-key-grad"><stop offset="0" stop-color="rgb(59,82,139)"/><stop offset="0.5" stop-color="rgb(253,231,37)"/><stop offset="1" stop-color="rgb(178,24,24)"/></linearGradient></defs><circle cx="8" cy="8" r="7" fill="url(#ap-key-grad)"/></svg>"""
const arrow* = """<svg width="26" height="10" viewBox="0 0 26 10" aria-hidden="true"><line x1="1" y1="5" x2="21" y2="5" style="stroke:var(--ap-ink)" stroke-width="1.5"/><path d="M19,1 L25,5 L19,9 z" style="fill:var(--ap-ink)"/></svg>"""
func legendHtml*(params: ApParameters): string =
## What color and arrows mean, then the parameters of the run.
let pref = if params.quantile == 0.5: "median similarity"
else: "similarity quantile " & $params.quantile
"<div class=\"key\">" &
"<span>" & evidenceDot & "</span><span>evidence r(k,k) + a(k,k) that k is an " &
"exemplar (it is one when positive)</span>" &
"<span>" & arrow & "</span><span>i → k: i's belief that k is its exemplar</span>" &
"</div>\n" &
"<p class=\"params\"><strong>parameters</strong> <span>preference " & pref & "</span> · <span>damping λ = " &
$params.damping & "</span> · <span>convits " & $params.convits & "</span></p>"
func figureHtml*(id: string, params: ApParameters): string =
## The frame an animation mounts in (the div `id`), with its legend.
"<figure>\n<div id=\"" & id & "\"></div>\n<figcaption>\n" & legendHtml(params) &
"\n</figcaption>\n</figure>"
when defined(js):
proc mountAnimation*(rootId: cstring, run: ApRun, px = 400.0, msPerStep = 90) =
let pts = run.input.points
let n = pts.len
let fr = frameOf(pts, px)
let last = run.iterations.r.len # 0 is before any message, then 1 .. last
let temperature = abs(run.preference) / 10
let zero = zeros[float](n)
let rad = 0.02 * px
# the evidence scale is the range of the whole run, so a color means the
# same at every iteration
var lo = Inf
var hi = -Inf
for t in 1 .. last:
for k in 0 ..< n:
let e = evidence(run.iterations.r[t - 1], run.iterations.a[t - 1], k)
lo = min(lo, e)
hi = max(hi, e)
# keep 0 strictly inside the scale, even for a run with no exemplar
lo = min(lo, -1e-9)
hi = max(hi, 1e-9)
# for how many iterations, up to t, each point has been an exemplar: points
# become opaque as they settle, over `convits` iterations
var streak = newSeq[seq[int]](last + 1)
streak[0] = newSeq[int](n)
for t in 1 .. last:
streak[t] = streak[t - 1]
for k in 0 ..< n:
if k in run.iterations.exemplars[t - 1]: inc streak[t][k]
else: streak[t][k] = 0
let settle = max(run.params.convits, 1)
let barH = 0.13 * px # the color bar under the points
var t = 0
var playing = false
var timer: Interval
var kxi: KaraxInstance # our own instance: several animations can share a page
proc stop() =
if playing:
clearInterval(timer)
playing = false
proc go(j: int) =
stop()
t = clamp(j, 0, last)
proc tick() =
if t < last: inc t
else: stop()
redraw(kxi)
proc toggle() =
if playing: stop()
else:
if t >= last: t = 0
playing = true
timer = setInterval(tick, msPerStep)
proc picture(): VNode =
let r = if t == 0: zero else: run.iterations.r[t - 1]
let a = if t == 0: zero else: run.iterations.a[t - 1]
let exemplars = if t == 0: newSeq[int]() else: run.iterations.exemplars[t - 1]
result = buildHtml(svg(viewBox = "0 0 " & f(px) & " " & f(px), class = "ap-plot",
role = "img")):
defs:
marker(id = "ap-head", viewBox = "0 0 10 10", refX = "9", refY = "5",
markerWidth = "5", markerHeight = "5", orient = "auto-start-reverse"):
path(d = "M0,0 L10,5 L0,10 z", class = "ap-head")
rect(width = f(px), height = f(px), class = "ap-frame")
for i in 0 ..< n:
let b = beliefs(r, a, i, temperature)
for k in 0 ..< n:
if k != i and b[k] >= 0.02:
let (x1, y1) = (sx(fr, pts[i].x), sy(fr, pts[i].y))
let (x2, y2) = (sx(fr, pts[k].x), sy(fr, pts[k].y))
let d = max(hypot(x2 - x1, y2 - y1), 1e-9)
let stopAt = rad * (if k in exemplars: 1.9 else: 1.4)
line(x1 = f(x1), y1 = f(y1), x2 = f(x2 - (x2 - x1) / d * stopAt),
y2 = f(y2 - (y2 - y1) / d * stopAt), class = "ap-arrow",
opacity = f(b[k]), `stroke-width` = f(0.6 + 1.4 * b[k]),
`marker-end` = "url(#ap-head)")
for k in 0 ..< n:
let e = evidence(r, a, k)
let isEx = k in exemplars
# before any message there is no evidence yet: a neutral point.
# Points are see-through until they become exemplars, then turn
# opaque as they stay exemplars (settled after `convits` iterations)
let color = if t == 0: "var(--muted)" else: evidenceColor(e, lo, hi)
let opacity = 0.35 + 0.65 * min(streak[t][k] / settle, 1.0)
circle(cx = f(sx(fr, pts[k].x)), cy = f(sy(fr, pts[k].y)),
r = f(if isEx: 1.5 * rad else: rad),
class = (if isEx: "ap-point ap-exemplar" else: "ap-point"),
fill = color, stroke = color, `fill-opacity` = f(opacity)):
title: text "point " & $k & ": evidence r(k,k) + a(k,k) = " &
formatFloat(e, ffDecimal, 3)
proc colorBar(): VNode =
## The evidence scale, linear from the lowest to the highest of the run,
## with 0 marked and a tick for every point now. The bar stretches to the
## figure's width; its labels are HTML, so they keep their size.
let r = if t == 0: zero else: run.iterations.r[t - 1]
let a = if t == 0: zero else: run.iterations.a[t - 1]
const (w, h, steps) = (1000.0, 24.0, 60)
proc bx(e: float): float = w * (e - lo) / (hi - lo)
result = buildHtml(tdiv(class = "ap-scale")):
tdiv(class = "ap-scale-row"):
span: text "← not an exemplar"
span: text "exemplar →"
svg(viewBox = "0 0 " & f(w) & " " & f(h), preserveAspectRatio = "none",
class = "ap-bar", role = "img", `aria-label` = "Evidence color scale"):
for j in 0 ..< steps:
# one rect per step, colored at its middle
let e = lo + (hi - lo) * (j.float + 0.5) / steps.float
rect(x = f(j.float * w / steps.float), y = "6", width = f(w / steps.float + 1),
height = "12", fill = evidenceColor(e, lo, hi))
for x in [0.5, bx(0), w - 0.5]:
line(x1 = f(x), y1 = "0", x2 = f(x), y2 = f(h), class = "ap-mark")
if t > 0:
for k in 0 ..< n:
let x = bx(evidence(r, a, k))
line(x1 = f(x), y1 = "2", x2 = f(x), y2 = "22", class = "ap-tick")
tdiv(class = "ap-scale-zero"):
span(style = style(StyleAttr.left, f(100 * bx(0) / w) & "%")): text "0"
proc createDom(): VNode =
let exemplars = if t == 0: newSeq[int]() else: run.iterations.exemplars[t - 1]
let label = if t == 0: "Before any message"
elif t == last and run.result.converged: "Iteration " & $t & ", converged"
else: "Iteration " & $t
let playLabel = if playing: "Pause" else: "Play"
result = buildHtml(tdiv(class = "ap-anim")):
picture()
colorBar()
tdiv(class = "ap-controls"):
button(class = "ap-play", onclick = toggle): text playLabel
button(class = "ap-step", title = "Previous iteration",
onclick = proc() = go(t - 1)): text "‹"
button(class = "ap-step", title = "Next iteration",
onclick = proc() = go(t + 1)): text "›"
input(`type` = "range", id = "ap-slider", min = "0", max = $last,
value = $t, class = "ap-slider", `aria-label` = "Iteration"):
proc oninput(ev: Event; n: VNode) = go(parseInt($n.value))
tdiv(class = "ap-readout"):
span(class = "ap-it"): text label
span(class = "ap-ex"):
text $exemplars.len & (if exemplars.len == 1: " exemplar" else: " exemplars")
kxi = setRenderer(createDom, root = rootId)
redraw(kxi)
proc enableArrowKeys*() =
## Left and right arrow keys step the page's first animation, except while
## typing in an input (the slider handles its own arrows).
document.addEventListener("keydown", proc (ev: Event) =
let e = KeyboardEvent(ev)
if ev.target != nil and Element(ev.target).nodeName == "INPUT": return
let steps = document.querySelectorAll(".ap-step")
if steps.len < 2: return
if e.key == "ArrowLeft":
steps[0].click()
ev.preventDefault()
elif e.key == "ArrowRight":
steps[1].click()
ev.preventDefault())
The page. The theme plugs our stylesheet into nimib, with no highlight.js, and adds
nbHeader, the title with the line above it and a link to the repository; the
stylesheet has the colors for light and dark.
theme.nim · 45 lines
## The look of the explanation pages: nimib's default theme with our own
## stylesheet (style.css) inlined, no highlight.js (code is highlighted at build
## time), so nothing loads from a CDN but the fonts; and nbHeader, the page's
## title block.
import std/json
import nimib
const fonts = """<link rel="preconnect" href="https://fonts.googleapis.com">
<link rel="stylesheet" href="https://fonts.googleapis.com/css2?family=Source+Serif+4:opsz,wght@8..60,400;8..60,600&family=Source+Sans+3:wght@400;600&family=JetBrains+Mono:wght@400&display=swap">"""
const style = "<style>\n" & staticRead("style.css") & "</style>"
func noHeader(blk: JsonNode, nb: Nb): string = ""
const githubLogo = """<svg aria-hidden="true" width="1.25em" height="1.25em" viewBox="0 0 16 16"><path fill="currentColor" fill-rule="evenodd" d="M8 0C3.58 0 0 3.58 0 8c0 3.54 2.29 6.53 5.47 7.59c.4.07.55-.17.55-.38c0-.19-.01-.82-.01-1.49c-2.01.37-2.53-.49-2.69-.94c-.09-.23-.48-.94-.82-1.13c-.28-.15-.68-.52-.01-.53c.63-.01 1.08.58 1.23.82c.72 1.21 1.87.87 2.33.66c.07-.52.28-.87.51-1.07c-1.78-.2-3.64-.89-3.64-3.95c0-.87.31-1.59.82-2.15c-.08-.2-.36-1.02.08-2.12c0 0 .67-.21 2.2.82c.64-.18 1.32-.27 2-.27c.68 0 1.36.09 2 .27c1.53-1.04 2.2-.82 2.2-.82c.44 1.1.16 1.92.08 2.12c.51.56.82 1.27.82 2.15c0 3.07-1.87 3.75-3.65 3.95c.29.25.54.73.54 1.48c0 1.07-.01 1.93-.01 2.2c0 .21.15.46.55.38A8.013 8.013 0 0 0 16 8c0-4.42-3.58-8-8-8z"></path></svg>"""
# The page's title, with an eyebrow above it (the line that says what the page
# belongs to) and, above it on the right, a GitHub link to the repository; it
# replaces nimib's header bar.
newNbBlock(NbHeader):
title: string
eyebrow: string
repo: string
toHtml:
"<header class=\"nb-header\">" &
(if blk.repo.len > 0: "<p class=\"repo\"><a href=\"https://github.com/" & blk.repo &
"\" title=\"source on GitHub\">" & githubLogo & "<span>" & blk.repo & "</span></a></p>" else: "") &
"<p class=\"eyebrow\">" & blk.eyebrow & "</p>" &
"<h1>" & blk.title & "</h1></header>"
proc header*(nb: var Nb, title: string, eyebrow = "", repo = "") =
nb.add newNbHeader(title = title, eyebrow = eyebrow, repo = repo)
nb.title = title
template nbHeader*(title: string, eyebrow = "", repo = "") = nb.header(title, eyebrow, repo)
proc useExplainTheme*(nb: var Nb) =
## Call right after `nbInit`.
nb.doc.context["stylesheet"] = %(fonts & "\n" & style)
nb.doc.context["highlight"] = %""
nb.doc.context["nb_style"] = %""
nb.doc.context["disableHighlightJs"] = %true
nb.backend.partials["header"] = noHeader
style.css · 89 lines
/* The stylesheet of the explanation pages, inlined by theme.nim: colors as
custom properties for light and dark, the page layout, code highlighting
(server side, highlight.js class names) and the animation's figure. */
:root {
--bg: #f6f7f9; --surface: #ffffff; --fg: #1c2330; --muted: #5d6675; --rule: #d9dee6;
--link: #2f6db0; --code-bg: #eef1f5;
--ap-ink: #1c2330; --ap-point: #2f6db0; --ap-frame: #ffffff;
--hl-keyword: #a626a4; --hl-string: #50a14f; --hl-number: #986801; --hl-comment: #6a737d;
--hl-type: #c18401; --hl-builtin: #0f7f8f;
--display: "Source Serif 4", Georgia, serif;
--body: "Source Sans 3", "Segoe UI", system-ui, sans-serif;
--mono: "JetBrains Mono", ui-monospace, Menlo, monospace;
}
@media (prefers-color-scheme: dark) { :root:not([data-theme="light"]) {
--bg: #11151b; --surface: #171c24; --fg: #e3e8ef; --muted: #98a2b1; --rule: #2a313c;
--link: #79a9e0; --code-bg: #1b212b;
--ap-ink: #e3e8ef; --ap-point: #79a9e0; --ap-frame: #171c24;
--hl-keyword: #c792ea; --hl-string: #a5d6a7; --hl-number: #f2b179; --hl-comment: #8b95a5;
--hl-type: #ffcb6b; --hl-builtin: #6fd3e0; color-scheme: dark } }
:root[data-theme="dark"] {
--bg: #11151b; --surface: #171c24; --fg: #e3e8ef; --muted: #98a2b1; --rule: #2a313c;
--link: #79a9e0; --code-bg: #1b212b;
--ap-ink: #e3e8ef; --ap-point: #79a9e0; --ap-frame: #171c24;
--hl-keyword: #c792ea; --hl-string: #a5d6a7; --hl-number: #f2b179; --hl-comment: #8b95a5;
--hl-type: #ffcb6b; --hl-builtin: #6fd3e0; color-scheme: dark }
body { background: var(--bg); color: var(--fg); font: 16px/1.55 var(--body); margin: 0; padding: 0 16px; }
main, footer { max-width: 44rem; margin: 0 auto; }
.nb-header { display: grid; gap: 8px; }
.nb-header .repo { margin: 0; max-width: none; text-align: right; font: 0.8rem/1 var(--mono); }
.nb-header .repo a { display: inline-flex; align-items: center; gap: 8px; color: var(--muted); text-decoration: none; white-space: nowrap; }
.nb-header .repo a:hover { color: var(--fg); }
.nb-header .repo svg { display: block; width: 1.25em; height: 1.25em; flex: none; }
main { padding-block: 32px; display: grid; gap: 16px; min-width: 0; }
main > * { min-width: 0; }
.eyebrow { font: 600 0.75rem/1 var(--body); letter-spacing: 0.08em; text-transform: uppercase; color: var(--muted); margin: 0; }
h1, h2 { font-family: var(--display); font-weight: 600; line-height: 1.2; margin: 0; text-wrap: balance; }
h1 { font-size: 1.75rem; }
h2 { font-size: 1.3rem; margin-top: 16px; }
p { margin: 0; max-width: 65ch; }
a { color: var(--link); }
code { font-family: var(--mono); font-size: 0.88em; }
pre { margin: 0; background: var(--code-bg); border: 1px solid var(--rule); border-radius: 6px; padding: 12px 14px; overflow-x: auto; font-size: 0.82rem; line-height: 1.5; }
.hljs-keyword, .hljs-meta { color: var(--hl-keyword); }
.hljs-string { color: var(--hl-string); }
.hljs-number, .hljs-literal { color: var(--hl-number); }
.hljs-comment { color: var(--hl-comment); font-style: italic; }
.hljs-type { color: var(--hl-type); }
.hljs-built_in { color: var(--hl-builtin); }
figure { margin: 0; background: var(--surface); border: 1px solid var(--rule); border-radius: 6px; padding: 16px; display: grid; gap: 12px; }
.ap-anim { display: grid; gap: 12px; }
.ap-plot { width: 100%; max-width: 100%; aspect-ratio: 1; display: block; }
.ap-frame { fill: var(--ap-frame); }
.ap-arrow { stroke: var(--ap-ink); }
.ap-head { fill: var(--ap-ink); }
.ap-point { stroke-width: 1.2; }
.ap-scale { display: grid; gap: 2px; }
.ap-scale-row { display: flex; justify-content: space-between; gap: 8px; font: 0.8rem var(--mono); color: var(--muted); }
.ap-scale-zero { position: relative; height: 1.2em; font: 0.8rem var(--mono); color: var(--muted); }
.ap-scale-zero span { position: absolute; transform: translateX(-50%); }
.ap-bar { width: 100%; height: 24px; display: block; }
.ap-tick, .ap-mark { vector-effect: non-scaling-stroke; }
.ap-tick { stroke: var(--ap-ink); stroke-width: 1.5; opacity: 0.75; }
.ap-mark { stroke: var(--muted); stroke-width: 1; }
.ap-controls { display: flex; align-items: center; gap: 8px; flex-wrap: wrap; }
button { font: 600 0.95rem var(--body); color: var(--fg); background: var(--bg); border: 1px solid var(--rule); border-radius: 4px; padding: 6px 12px; cursor: pointer; }
button:hover { border-color: var(--muted); }
button:focus-visible, .ap-slider:focus-visible { outline: 2px solid var(--ap-point); outline-offset: 2px; }
.ap-play { min-width: 5.5em; }
.ap-step { font-size: 1.1rem; padding: 4px 12px; }
.ap-slider { flex: 1 1 12rem; min-width: 0; accent-color: var(--ap-point); }
.ap-readout { display: flex; justify-content: space-between; gap: 12px; flex-wrap: wrap; font: 0.9rem var(--mono); font-variant-numeric: tabular-nums; color: var(--muted); }
figcaption { font-size: 0.95rem; color: var(--muted); display: grid; gap: 8px; border-top: 1px solid var(--rule); padding-top: 12px; }
.key { display: grid; grid-template-columns: 28px 1fr; gap: 6px 10px; align-items: start; }
.key > span:nth-child(odd) { display: flex; justify-content: center; height: 1.5em; align-items: center; }
.params { max-width: none; font: 0.85rem var(--mono); border-top: 1px solid var(--rule); padding-top: 8px; }
.params span { white-space: nowrap; }
.params strong { font: 600 0.85rem var(--body); margin-right: 4px; }
details.nb-file { border: 1px solid var(--rule); border-radius: 6px; background: var(--surface); }
details.nb-file > summary { cursor: pointer; padding: 10px 14px; font: 0.9rem var(--mono); }
details.nb-file > summary span { color: var(--muted); }
details.nb-file[open] > summary { border-bottom: 1px solid var(--rule); }
details.nb-file pre { border: 0; border-radius: 0 0 6px 6px; }
footer { border-top: 1px solid var(--rule); padding-block: 16px 32px; display: grid; gap: 16px; }
.nb-box { display: flex; align-items: center; justify-content: space-between; gap: 12px; flex-wrap: wrap; }
.nb-small { font-size: 0.85rem; color: var(--muted); }
button.nb-small { font-size: 0.85rem; color: var(--fg); padding: 4px 10px; }
section#source { display: none; }
@media (prefers-reduced-motion: reduce) { * { transition: none !important; } }