SPINE — User Guide¶
Scale Partitioning by Incremental Nested Entropy.
SPINE takes a nonnegative measure over a collection of objects and cuts it into scale strata. It finds where the collection separates into distinct scales, and how strongly, without being told how many groups to look for.
You supply two things: a measure, one nonnegative number per object, and an orientation, saying which end of that measure is the atypical one. Those are the scientific input. Everything after them is computed.
Run this guide cell by cell. Sections 1 to 5 are the basics. Sections 6 onward are worked scenarios, each on a different kind of data, and they can be read in any order.
1. Before you start¶
pip install entropy-spine — or, inside the repository, add python/src to
your path. Nothing else to set up.
%matplotlib inline
import numpy as np
from entropy_spine import (spine, spine_summary, spine_compare,
spine_data_plot, spine_plot, spine_network_plot,
spine_version)
spine_version();
entropy-spine 0.1.0 package .../python/src/entropy_spine python 3.13.14 numpy 2.5.0 path clean, no shadowing copies
2. Your first call¶
Two hundred objects drawn from a standard normal cloud, plus twelve sitting a long way out. Nothing tells SPINE that there are two groups.
rng = np.random.default_rng(11)
X = np.vstack([rng.standard_normal((200, 2)),
rng.standard_normal((12, 2)) + 9])
scores, result = spine(X, "dist2mean", "upper")
print(f"strata found {result.n_strata}")
print(f"cut at measure {result.boundary_values[0]:.4g}")
print(f"objects flagged {int(np.sum(scores > 0))} of {scores.size}")
strata found 2 cut at measure 3.423 objects flagged 12 of 212
3. Reading the result¶
Scores run from 0 for the typical stratum to 1 for the most atypical, one
per object, in your object order. result carries the rest, also in your
order, and spine_summary says it in words. It recomputes nothing, so it
also works on a Result loaded back from a file long after the run.
spine_summary(result);
SPINE 212 objects, 2 strata measure dist2mean | orientation upper
stratum objects measure range score
1 200 0.01647 to 3.423 0.000
2 12 10.96 to 13.68 1.000
cut at 3.423
result.boundary_values applies these cuts to new data.
result.boundary_values is the result you can carry elsewhere. It is a value
in your own measure's units, so the same cut can be applied to new data or
quoted in a paper. result.info["boundaries"], by contrast, are ranks into
an internally sorted copy and mean nothing outside this call.
4. Seeing it¶
Two figures, each with its own job.
spine_data_plot draws your data in the space it lives in, coloured by
score. Objects are drawn in increasing order of score, so the atypical ones
are never buried under the bulk, and the colour bar is ticked at the strata
with their populations, so it doubles as the legend.
spine_data_plot(result, title="a bulk and a far group");
spine_plot shows where the cuts came from: the sorted measure with the cut marked, and the incremental nested entropy profile behind it.
spine_plot(result, title="a bulk and a far group");
5. The measure is the modelling step¶
This is the part worth slowing down for. The same data under two measures can give opposite answers, and both are correct — they answer different questions.
Three groups: a bulk, a halo displaced a little and spread widely, and an island displaced a long way and packed tight.
rng = np.random.default_rng(4242)
bulk = rng.standard_normal((200, 2))
ang = 2 * np.pi * rng.random(30)
halo = 5.0 * np.column_stack([np.cos(ang), np.sin(ang)]) \
+ 2.0 * rng.standard_normal((30, 2))
island = [12, 12] + 0.25 * rng.standard_normal((20, 2))
D = np.vstack([bulk, halo, island])
spine_compare(D, ["dist2mean", "robustrownorm", "knn", "dist2median"],
"upper", k=5);
SPINE comparison over 4 measures
measure strata flagged fraction first cut
dist2mean 2 27 0.108 6.485
robustrownorm 2 42 0.168 3.157
knn 3 32 0.128 1.071
dist2median 2 41 0.164 4.157
agreement on which objects are atypical (Jaccard)
dist2mean robustrownorm knn dist2median
dist2mean - 0.64 0.13 0.66
robustrownorm 0.64 - 0.42 0.98
knn 0.13 0.42 - 0.40
dist2median 0.66 0.98 0.40 -
dist2mean and knn overlap on only 13% of what they flag.
They are answering different questions. Deciding which one you
meant is the modelling step, not a detail. See the measure guide.
Read the island through the agreement matrix. dist2mean calls it maximally
atypical. knn calls it perfectly typical and picks out the halo instead:
the island is far away but densely packed, so its nearest-neighbour
distances are small.
Distance from a centre and local density are different questions, and the measure is where you say which one you meant. Neither answer is wrong. The measure guide on the documentation site is the fuller treatment.
Side by side, the difference is visible rather than tabulated.
_, cmp_detail = spine_compare(D, ["dist2mean", "knn"], "upper", k=5,
report=False)
spine_data_plot(cmp_detail["results"][0],
title="dist2mean: far from the centre");
spine_data_plot(cmp_detail["results"][1],
title="knn: in a sparse neighbourhood");
6. Scenario: rare categories¶
Data structure: a plain vector of counts.
Forty tokens seen once or twice, six seen dozens of times. Here the RARE
ones are the interesting ones, so small is atypical and the orientation is
"lower".
rng = np.random.default_rng(20)
counts = np.concatenate([rng.integers(1, 4, 40),
rng.integers(55, 71, 6)]).astype(float)
s_freq, res_freq = spine(counts, "frequency", "lower")
print(f"strata {res_freq.n_strata}, flagged "
f"{int(np.sum(s_freq > 0))} of {s_freq.size} as rare")
print(f"cut at frequency {res_freq.boundary_values[0]:.4g}")
strata 2, flagged 40 of 46 as rare cut at frequency 0.006652
Getting the orientation backwards does not produce an error. It produces a confident answer to a question you did not mean to ask.
6b. Scenario: rare states in symbolic data¶
Data structure: text — a sequence of strings, or a DataFrame mixing text and numbers.
Use frequency when you already hold one count per category. Use
statefrequency when you hold raw observations that still need tallying:
each ROW is one observed state, rows are matched exactly, and the objects
become the distinct states.
Four thousand draws from a sixty-word vocabulary with a heavy head.
rng = np.random.default_rng(9)
n_tok = 60
weights = 1.0 / (np.arange(1, n_tok + 1) ** 1.1)
draws = rng.choice(n_tok, 4000, p=weights / weights.sum())
words = np.asarray([f"w{d:02d}" for d in draws])
s_word, res_word = spine(words, "statefrequency", "lower")
print(f"{draws.size} observations over {s_word.size} distinct states")
print(f"strata {res_word.n_strata}, states flagged as rare "
f"{int(np.sum(s_word > 0))}")
4000 observations over 60 distinct states strata 2, states flagged as rare 49
Text and numbers can be mixed, and the state is the whole row.
import pandas as pd
log = pd.DataFrame({"status": ["ok", "ok", "warn", "ok",
"fault", "warn", "ok", "fault"],
"code": [0, 0, 1, 0, 2, 1, 0, 3]})
_, res_log = spine(log, "statefrequency", "lower")
list(zip(res_log.object_ids, res_log.measure["counts"].tolist()))
[(np.str_('fault, 2'), 1),
(np.str_('fault, 3'), 1),
(np.str_('ok, 0'), 4),
(np.str_('warn, 1'), 2)]
Two cautions. Numeric columns are matched at full precision, so round or bin first if near-misses should count as the same state. And the object count is the number of DISTINCT states, not of observations: forty thousand rows over three states give SPINE three objects, which is not enough to partition.
To score the original rows rather than the states:
row_scores = s_word[res_word.measure["stateOfRow"] - 1]
7. Scenario: hubs in a network¶
Data structure: a square adjacency matrix.
The objects are nodes, the measure is the degree, and the orientation is
"upper" because hubs are the atypical ones. Three networks, each asking a
different question of the same measure. The generators are defined at the
end of this section.
A. A small-world backbone with five designated hubs. Tight, regular degrees plus a handful of deliberate hubs, so the ground truth is unambiguous.
def small_world(n, k, beta, seed):
"""Watts-Strogatz: a ring lattice with each edge rewired with
probability beta. Degrees stay tightly concentrated around k."""
rng = np.random.default_rng(seed)
A = np.zeros((n, n))
for i in range(n):
for j in range(1, k // 2 + 1):
t = (i + j) % n
A[i, t] = A[t, i] = 1
for i in range(n):
for j in range(1, k // 2 + 1):
t = (i + j) % n
if A[i, t] == 1 and rng.random() < beta:
candidates = np.flatnonzero(A[i] == 0)
candidates = candidates[candidates != i]
if candidates.size:
A[i, t] = A[t, i] = 0
nt = rng.choice(candidates)
A[i, nt] = A[nt, i] = 1
return A
def add_hubs(A, n_hub, frac, seed):
"""Attach n_hub new nodes, each wired to a fraction of the network."""
rng = np.random.default_rng(seed)
n_base = A.shape[0]
out = np.zeros((n_base + n_hub, n_base + n_hub))
out[:n_base, :n_base] = A
hub_ids = np.arange(n_base, n_base + n_hub)
for v in hub_ids:
targets = rng.permutation(n_base)[: round(frac * n_base)]
out[v, targets] = out[targets, v] = 1
return out, hub_ids
W = small_world(200, 6, 0.10, 7)
A1, hub_ids1 = add_hubs(W, 5, 0.18, 7)
truth1 = np.zeros(A1.shape[0], dtype=bool)
truth1[hub_ids1] = True
s1, r1 = spine(A1, "degree", "upper")
deg1 = A1.sum(axis=1)
print(f"A. {A1.shape[0]} nodes, degrees {int(deg1.min())} to "
f"{int(deg1.max())}, cut at {r1.boundary_values[0]:g}")
print(f" designated hubs found {int(np.sum((s1 > 0) & truth1))} of "
f"{int(truth1.sum())}, other nodes flagged "
f"{int(np.sum((s1 > 0) & ~truth1))}")
A. 205 nodes, degrees 5 to 36, cut at 11 designated hubs found 5 of 5, other nodes flagged 0
B. A scale-free network with six designated hubs. Preferential attachment makes hubs of its own accord, so the designated ones are not the only real hubs present.
def scale_free(n, m, seed):
"""Barabasi-Albert preferential attachment."""
rng = np.random.default_rng(seed)
m0 = m + 1
A = np.zeros((n, n))
A[:m0, :m0] = 1 - np.eye(m0)
pool = list(np.repeat(np.arange(m0), m0 - 1))
for v in range(m0, n):
chosen = []
while len(chosen) < m:
t = pool[rng.integers(len(pool))]
if t not in chosen:
chosen.append(t)
for t in chosen:
A[v, t] = A[t, v] = 1
pool.extend(chosen)
pool.extend([v] * m)
return A
S = scale_free(180, 2, 20260810)
A2, hub_ids2 = add_hubs(S, 6, 0.22, 20260810)
truth2 = np.zeros(A2.shape[0], dtype=bool)
truth2[hub_ids2] = True
s2, r2 = spine(A2, "degree", "upper")
deg2 = A2.sum(axis=1)
extra = np.flatnonzero((s2 > 0) & ~truth2)
print(f"B. designated hubs found {int(np.sum((s2 > 0) & truth2))} of "
f"{int(truth2.sum())}, other nodes flagged {extra.size}")
print(f" those extra nodes have degrees "
f"{sorted(int(d) for d in deg2[extra])}, against a bulk maximum of "
f"{int(deg2[s2 == 0].max())}")
B. designated hubs found 6 of 6, other nodes flagged 6 those extra nodes have degrees [14, 20, 26, 31, 35, 53], against a bulk maximum of 11
They are not false positives. Preferential attachment produced genuine hubs, and the degree measure does not care which were designed and which grew.
C. The same backbone with no hubs at all. One stratum is the answer: these degrees sit at a single scale. That is a property of the network, not of the cut.
A3 = small_world(200, 6, 0.10, 7)
s3, r3 = spine(A3, "degree", "upper")
deg3 = A3.sum(axis=1)
print(f"C. degrees {int(deg3.min())} to {int(deg3.max())}, "
f"strata {r3.n_strata}, flagged {int(np.sum(s3 > 0))}")
C. degrees 4 to 8, strata 1, flagged 0
spine_network_plot draws the graph itself: a force layout with nodes coloured by score and sized by degree, so hubs read as hubs.
spine_network_plot(r1, A1, title="small-world with designated hubs", seed=3);
spine_network_plot(r2, A2, title="scale-free", seed=3);
8. Scenario: how many components does this matrix have?¶
Data structure: a rectangular data matrix, treated spectrally.
Three strong directions buried in noise. The measure is the singular values themselves, so the OBJECTS are the singular values, not the rows.
rng = np.random.default_rng(77)
U = np.linalg.qr(rng.standard_normal((60, 60)))[0]
Vt = np.linalg.qr(rng.standard_normal((25, 25)))[0]
sv = np.concatenate([[40, 26, 15], 0.6 * np.ones(22)])
M = U[:, :25] @ np.diag(sv) @ Vt.T
s_sv, res_sv = spine(M, "singular", "upper")
print(f"objects (singular values) {s_sv.size}")
print(f"strata {res_sv.n_strata}, components above the noise floor: "
f"{int(np.sum(s_sv > 0))}")
objects (singular values) 25 strata 2, components above the noise floor: 3
Because the objects are no longer the rows, there is no per-object data to draw, and spine_data_plot would tell you so. spine_plot still works, and the measure panel is exactly the spectrum with the cut on it.
spine_plot(res_sv.info, title="singular values, with the cut");
9. Scenario: a DataFrame with identifiers¶
Data structure: a pandas DataFrame with a meaningful index.
Fifty-five ordinary sites and five anomalous ones. The identifiers survive the round trip, so you get back the names of what was flagged, not row numbers.
rng = np.random.default_rng(8)
sites = np.vstack([rng.standard_normal((55, 3)),
rng.standard_normal((5, 3)) + 7])
T = pd.DataFrame(sites, columns=["x", "y", "z"],
index=[f"site{i:02d}" for i in range(1, 61)])
s_t, res_t = spine(T, "robustrownorm", "upper")
print(f"table: {res_t.mu.size} objects, {res_t.n_strata} strata")
print([res_t.object_ids[i] for i in np.flatnonzero(s_t > 0)])
table: 60 objects, 2 strata ['site56', 'site57', 'site58', 'site59', 'site60']
10. Scenario: a sequence of records¶
Data structure: a list of dicts, one per object.
Numeric values become the features. A text value becomes the identifiers, automatically.
rng = np.random.default_rng(15)
parts = []
for i in range(1, 46):
far = i > 41
parts.append({
"label": f"part-{i:03d}",
"width": 10 + rng.standard_normal() + 14 * far,
"height": 5 + rng.standard_normal() + 9 * far,
"mass": 2 + rng.standard_normal() + 6 * far,
})
s_s, res_s = spine(parts, "robustrownorm", "upper")
print(f"records: {res_s.mu.size} objects, {res_s.n_strata} strata")
print([res_s.object_ids[i] for i in np.flatnonzero(s_s > 0)])
records: 45 objects, 2 strata ['part-042', 'part-043', 'part-044', 'part-045']
11. Scenario: three dimensions¶
Data structure: an n-by-3 matrix.
A dense core, a thin shell around it, and far outliers.
rng = np.random.default_rng(90210)
V1 = rng.standard_normal((90, 3))
U1 = V1 / np.linalg.norm(V1, axis=1, keepdims=True)
V2 = rng.standard_normal((18, 3))
U2 = V2 / np.linalg.norm(V2, axis=1, keepdims=True)
X3 = np.vstack([0.9 * rng.standard_normal((260, 3)),
U1 * (7 + 0.25 * rng.standard_normal((90, 1))),
U2 * (34 + 3 * rng.standard_normal((18, 1)))])
_, res3 = spine(X3, "dist2median", "upper")
cuts = ", ".join(f"{c:.3g}" for c in res3.boundary_values)
print(f"3-D: {res3.n_strata} strata, cut radii {cuts}")
3-D: 3 strata, cut radii 3.06, 7.56
The three strata score 0, 0.18 and 1, so on a score colour scale the core and the shell come out nearly the same shade. Colouring by stratum spaces them evenly instead, and the bar is labelled "stratum" so nothing is misread.
spine_data_plot(res3, title="core, shell, outliers", colour_by="stratum");
12. Scenario: many features at once¶
Data structure: a wide matrix, more features than can be drawn.
Twelve features. The measure reduces them to one number as usual; only the drawing is affected, and it shows the first two principal components with the axes labelled to say so.
rng = np.random.default_rng(64)
X12 = np.vstack([rng.standard_normal((240, 12)),
rng.standard_normal((18, 12)) + 5])
_, res12 = spine(X12, "robustrownorm", "upper")
print(f"12 features: {res12.n_strata} strata, "
f"{int(np.sum(res12.scores > 0))} flagged")
spine_data_plot(res12, title="12 features, projected");
12 features: 2 strata, 18 flagged
13. Scenario: you already have the measure¶
Data structure: a vector that is already the quantity you care about.
If your domain hands you the number — an energy, a claim size, a residual —
use it directly with asis. A domain quantity almost always beats a
constructed one, because it carries meaning your reader already understands.
rng = np.random.default_rng(5)
energy = np.concatenate([rng.gamma(2, 1.0, 300), rng.gamma(2, 25, 9)])
s_e, res_e = spine(energy, "asis", "upper")
print(f"asis: {res_e.n_strata} strata, cut at "
f"{res_e.boundary_values[0]:.4g}, {int(np.sum(s_e > 0))} flagged")
asis: 2 strata, cut at 6.299, 10 flagged
14. Scenario: atypical at both ends¶
Sometimes both extremes matter and the typical objects are in the middle. Bilateral orientation needs an anchor: an object you consider typical.
Here the measure is distance from the median, and the anchor is the object closest to it.
rng = np.random.default_rng(42)
v = np.concatenate([rng.standard_normal(200) + 20,
0.5 * rng.standard_normal(8) + 3,
0.5 * rng.standard_normal(8) + 46])
dev = np.abs(v - np.median(v))
anchor = int(np.argmin(dev)) + 1 # object ids default to 1..N
s_b, res_b = spine(dev, "asis", "bilateral", bilateral_anchor=anchor)
print(f"bilateral: {res_b.n_strata} strata, "
f"{int(np.sum(s_b > 0))} flagged at the two ends")
print(f"anchor object {anchor} has score {s_b[anchor - 1]:.3f}")
bilateral: 2 strata, 17 flagged at the two ends anchor object 20 has score 0.000
15. Scenario: robust or not¶
When you are hunting for objects at a different scale, do not standardize by statistics those same objects distort. The standard deviation is inflated by the very group you are looking for; the MAD is not.
rng = np.random.default_rng(20260810)
Xr = np.vstack([rng.standard_normal((240, 8)),
rng.standard_normal((45, 8)) + 6,
rng.standard_normal((15, 8)) + 30])
_, rz = spine(Xr, "zrownorm", "upper")
_, rmad = spine(Xr, "robustrownorm", "upper")
print(f"zrownorm {rz.n_strata} strata")
print(f"robustrownorm {rmad.n_strata} strata (the planted structure)")
zrownorm 2 strata robustrownorm 3 strata (the planted structure)
Prefer robustrownorm over zrownorm, and dist2median over dist2mean, unless you have a reason not to.
16. Scenario: nothing to find¶
Structureless data gives one stratum, and that is an answer rather than a failure: this measure does not separate these objects into scales. It is a property of the measure you chose, not of the cut.
rng = np.random.default_rng(3)
_, res_flat = spine(rng.standard_normal((300, 4)), "dist2median", "upper")
print(f"unstructured data: {res_flat.n_strata} stratum")
unstructured data: 1 stratum
The right response is to ask whether a different measure would express the structure better, not to loosen anything — there is nothing to loosen.
17. Options you may want¶
objects="columns" your matrix has objects in columns, not rows
k=8 neighbour count for the knn measure
object_ids=ids supply identifiers yourself
missing="omit" drop objects with missing values; scores return NaN
thin=True very large collections: cut located on a thinned copy
boundary_mode="selected" keep a selected subset of the candidate cuts
bilateral_anchor=id the typical object, for bilateral orientation
profile_weights=w weight the profile ranks, selected mode onlyAn example of boundary_mode: candidate mode keeps every candidate
boundary, selected mode chooses among them.
rng = np.random.default_rng(99)
Xs = np.vstack([rng.standard_normal((150, 2)),
rng.standard_normal((30, 2)) + 4,
rng.standard_normal((10, 2)) + 14])
_, rc = spine(Xs, "dist2median", "upper")
_, rs = spine(Xs, "dist2median", "upper", boundary_mode="selected")
print(f"candidate mode {rc.n_strata} strata, "
f"selected mode {rs.n_strata} strata")
candidate mode 2 strata, selected mode 1 strata
18. Where to go next¶
- the measure guide on the documentation site — how to choose a measure, in full
- the repository README — the reference for everything
Resultgives you paper/figures/python/— the manuscript figures, reproduced in Python
help(spine_measure) lists the ten measures. help(spine) documents every
option. help(spine_data_plot), help(spine_network_plot) and
help(spine_plot) cover the figures.