VariableMarkovOracle.kt
src/jvmMain/kotlin/afield/math/oracle/VariableMarkovOracle.kt · 297 lines
package afield.math.oracle
import kotlin.math.ln
import kotlin.math.max
import kotlin.math.minBuilding a Variable Markov Oracle over a recording, and asking it where a performance could continue from.
VariableMarkovOracle indexes a Trajectory once, in its constructor, and then answers three
kinds of question about it. Its suffix links (sfx) and repeated-suffix lengths (lrs) say
where the material repeats itself, which is what Improviser walks. cue ranks positions by
how well the samples before them match a window of recent history, over any subset of channels,
which is how a performer's input finds a place in memory. bridges finds the nearest moments
in other episodes, where suffix links rarely reach, so several separately recorded takes can
be played as one memory.
The one parameter that matters is the threshold at which two frames count as the same symbol,
and build chooses it rather than asking for it: it constructs an oracle at each of a range of
candidate thresholds and keeps the one with the highest total information rate. That rate is
near zero at both ends of the range, since a threshold too small finds no repeats and one too
large makes everything a repeat, so the choice is a peak rather than a limit.
The construction is a Kotlin port of Cheng-i Wang's reference Python vmo package, kept in
Kotlin so an oracle can be built and queried inside the frame loop.
VmoReferenceTest checks the suffix links, the repeat lengths and the summed information rate
against a fixture generated from the pinned package.
class VariableMarkovOracle internal constructor
A Variable Markov Oracle (Wang and Dubnov, 2014) over a Trajectory: a Factor Oracle (Allauzen, Crochemore and Raffinot, 1999) built on continuous frames, where two frames count as "the same symbol" when their GestureFeatures lie within threshold of each other.
State i (1-based) holds sample i - 1; state 0 is the empty start. The structure is three
arrays:
- forward transitions — i → i + 1 along the recording, plus oracle transitions to later
states that continue the same context;
- suffix links sfx — the earlier state ending the longest repeated suffix of the material
up to i;
- lrs — the length of that repeated suffix, the amount of shared context.
Suffix links are what make recombination possible: jumping from i to sfx[i](../../../../../../reference/source/jvmMain/afield/audio/synth/Subharmonic.md#L240) and carrying
on from there continues something that, for the last lrs[i](../../../../../../reference/source/jvmMain/afield/audio/synth/Subharmonic.md#L240) samples, sounded the same.
Episode boundaries are respected: the first sample of an episode gets no forward transition from the last sample of the previous one, and lrs never reaches back across a boundary.
This is a Kotlin port of the construction in Cheng-i Wang's reference Python vmo package
(GPL-3.0, compatible with Field's AGPL-3.0+), kept in Kotlin so it can be stepped inside the
frame loop. VmoReferenceTest cross-checks the suffix links against the pinned package.
class VariableMarkovOracle internal constructor(
/** The remembered material. */
val trajectory: Trajectory,
/** The scaling and feature construction used for every comparison. */
val features: GestureFeatures,
/** Feature frames, one per sample. */
val frames: Array<DoubleArray>,
/** Frames closer than this are the same symbol. */
val threshold: Double,
) {
private val n = frames.size
private val trn = Array(n + 1) { IntArrayList() }
private val rsfx = Array(n + 1) { IntArrayList() }val sfx = IntArray(n + 1)
Suffix link per state; sfx[0] = -1.
val lrs = IntArray(n + 1)
Length of the repeated suffix per state.
val symbol = IntArray(n + 1)
Symbol (cluster of frames within threshold) per state; state 0 is symbol 0.
private var symbolCount = 1init
init {
sfx[0] = -1
for (i in 1..n) add(i)
}val size: Int get() = n + 1
Number of states (samples + 1).
private fun episodeStartState(i: Int): Int = trajectory.episodeStarts[trajectory.episodeOf(i - 1)] + 1
private fun add(i: Int) {
val x = frames[i - 1]
val newEpisode = i > 1 && trajectory.startsEpisode(i - 1)
var k: Int
var pi1: Int
if (newEpisode) { k = 0; pi1 = 0 } else { trn[i - 1].add(i); k = sfx[i - 1]; pi1 = i - 1 }
var best = -1
while (k != -1) {
var bd = Double.POSITIVE_INFINITY
trn[k].forEach { t -> if (t < i) { val d = distance(frames[t - 1], x); if (d < threshold && d < bd) { bd = d; best = t } } }
if (best >= 0) break
trn[k].add(i); pi1 = k; k = sfx[k]
}
if (k == -1 || best < 0) { sfx[i] = 0; lrs[i] = 0; symbol[i] = symbolCount++ }
else {
sfx[i] = best
lrs[i] = commonSuffix(pi1, best - 1) + 1
symbol[i] = symbol[best]
}
// the reference's "_find_better": prefer a sibling whose context extends one further
val want = symbol[i - lrs[i]]
val siblings = rsfx[sfx[i]]
for (q in 0 until siblings.size) {
val j = siblings[q]
if (lrs[j] == lrs[i] && j - lrs[i] >= 0 && symbol[j - lrs[i]] == want) { lrs[i]++; sfx[i] = j; break }
}
lrs[i] = min(lrs[i], i - episodeStartState(i) + 1)
rsfx[sfx[i]].add(i)
}
private fun commonSuffix(p1: Int, p2In: Int): Int {
var p2 = p2In
if (p2 == sfx[p1]) return lrs[p1]
while (p2 != 0 && sfx[p2] != sfx[p1]) p2 = sfx[p2]
return min(lrs[p1], lrs[p2])
}fun transitions(i: Int): IntArray = trn[i].toArray()
Forward transitions out of state i.
fun reverseSuffixes(i: Int): IntArray = rsfx[i].toArray()
States whose suffix link points at i.
val alphabetSize: Int get() = (1..n).count { sfx[it] == 0 }
Number of distinct symbols: states whose suffix link is the start state.
fun informationRate(): DoubleArray
Information rate per sample, in bits: the reference package's IR(ir_type='cum').
The oracle is compressed into Compror codewords (a new symbol, or a copy of a repeated
block). At sample i, h0 is log2 of the number of distinct symbols so far and h1 is
log2 of the number of codewords so far divided by the length of the block containing i.
The rate is max(0, h0 - h1): how much cheaper the past makes the present.
Too small a threshold makes every frame a new symbol (nothing repeats, rate ≈ 0); too
large makes everything one symbol (nothing to say, rate ≈ 0). build keeps the threshold
whose totalInformationRate is highest, as the reference's find_threshold does.
fun informationRate(): DoubleArray {
val cw0 = DoubleArray(n); val cw1 = DoubleArray(n); val bl = DoubleArray(n)
var j = 0; var i = 0; var at = 0
while (j < n) {
while (i < n && lrs[i + 1] >= i - j + 1) i++
if (i == j) { i++; cw0[at] = 1.0; cw1[at] = 1.0; bl[at] = 1.0; at++ }
else { val l = i - j; cw1[at] = 1.0; for (q in at until min(n, at + l)) bl[q] = l.toDouble(); at += l }
j = i
}
val out = DoubleArray(n)
var c0 = 0.0; var c1 = 0.0
val log2 = ln(2.0)
for (q in 0 until n) {
c0 += cw0[q]; c1 += cw1[q]
val h0 = if (c0 > 0) ln(c0) / log2 else Double.NEGATIVE_INFINITY
val h1 = if (c1 > 0) ln(c1) / log2 / bl[q] else Double.NEGATIVE_INFINITY
out[q] = max(0.0, h0 - h1).let { if (it.isNaN()) 0.0 else it }
}
return out
}fun totalInformationRate(): Double = informationRate().sum()
Sum of informationRate over the oracle: the quantity the threshold search maximises.
fun cue(history: List<DoubleArray>, window: Int = min(history.size, (trajectory.sampleRate).toInt()), mask: BooleanArray = BooleanArray(features.width) { true }
Ranked places in memory whose preceding window samples best match history, the
most recent samples last (historyt``, raw values).
Only channels whose mask entry is true are compared, so a cue can be made from the
channels a performer is driving while the rest are free. Candidates whose window would
cross an episode boundary are skipped. Distances are RMS per compared feature (scaled value or
velocity), so they do not grow with the number of channels. Each Cue's position is the sample to continue
from: the memory's next sample after the match is position + 1.
fun cue(history: List<DoubleArray>, window: Int = min(history.size, (trajectory.sampleRate).toInt()), mask: BooleanArray = BooleanArray(features.width) { true }, count: Int = 5, separation: Int = window / 2): List<Cue> {
val w = min(window, history.size)
require(w >= 1) { "cue needs at least one sample of history" }
val m = Matcher(history, w, mask)
val scored = ArrayList<Cue>()
for (end in w - 1 until n) {
val d = m.at(end)
if (d.isFinite()) scored.add(Cue(end, d, lrs[end + 1]))
}
scored.sortBy { it.distance }
val picked = ArrayList<Cue>()
for (c in scored) {
if (picked.size >= count) break
if (picked.none { kotlin.math.abs(it.position - c.position) < separation }) picked.add(c)
}
return picked
}
private val bridgeTable: Array<IntArray> by lazy {
val starts = trajectory.episodeStarts
if (starts.size < 2) return@lazy Array(n) { IntArray(0) }
val k = 3
Array(n) { i ->
val ep = trajectory.episodeOf(i)
val best = IntArray(k) { -1 }; val bd = DoubleArray(k) { Double.POSITIVE_INFINITY }
for (j in 0 until n) {
if (trajectory.episodeOf(j) == ep || trajectory.endsEpisode(j)) continue
val d = distance(frames[i], frames[j])
if (d < bd[k - 1]) {
var at = k - 1
while (at > 0 && d < bd[at - 1]) { bd[at] = bd[at - 1]; best[at] = best[at - 1]; at-- }
bd[at] = d; best[at] = j
}
}
best.filter { it >= 0 }.toIntArray()
}
}fun bridges(sample: Int): IntArray = bridgeTable[sample]
Cross-episode transitions from sample sample: up to three samples in other episodes whose
feature frames (value and velocity) are nearest to it, nearest first. The oracle's own suffix
links rarely cross between separately recorded takes; these give an improviser a way across.
Computed for every sample on first use (all pairs of samples, so a few seconds for a few
thousand samples). Empty for a single episode.
fun matchDistance(history: List<DoubleArray>, position: Int, window: Int = min(history.size, (trajectory.sampleRate).toInt()), mask: BooleanArray = BooleanArray
How well history matches the window samples of memory ending at position (the same
distance cue ranks by), or +∞ if that window crosses an episode boundary.
fun matchDistance(history: List<DoubleArray>, position: Int, window: Int = min(history.size, (trajectory.sampleRate).toInt()), mask: BooleanArray = BooleanArray(features.width) { true }): Double =
if (position < 0 || position >= n) Double.POSITIVE_INFINITY else Matcher(history, min(window, history.size), mask).at(position)
private inner class Matcher(history: List<DoubleArray>, val w: Int, mask: BooleanArray) {
private val h = history.subList(history.size - w, history.size)
private val hf = Array(w) { k -> features.frame(h[k], if (k > 0) h[k - 1] else null) }
private val dims = BooleanArray(features.dimension) { d -> mask[d % features.width] }
private val firstVelocity = features.dimension - features.width
fun at(end: Int): Double {
val start = end - w + 1
if (start < 0 || trajectory.episodeOf(start) != trajectory.episodeOf(end)) return Double.POSITIVE_INFINITY
var s = 0.0; var count = 0
for (k in 0 until w) {
val a = hf[k]; val b = frames[start + k]
// the first history sample has no velocity: compare its values only
for (d in a.indices) if (dims[d] && (k > 0 || d < firstVelocity)) { val e = a[d] - b[d]; s += e * e; count++ }
}
// RMS per compared feature, so distances mean the same for two channels or a hundred
return kotlin.math.sqrt(s / kotlin.math.max(1, count))
}
}override fun toString()
override fun toString() = "VariableMarkovOracle(${n} samples, threshold=${"%.3g".format(threshold)}, alphabet=$alphabetSize, IR=${"%.1f".format(totalInformationRate())} bits)"companion object
companion object {fun build(t: Trajectory, features: GestureFeatures = GestureFeatures.fit(t), threshold: Double? = null, candidates: Int = 24): VariableMarkovOracle
Builds the oracle for t. With no threshold it scans candidates thresholds and keeps
the one with the highest totalInformationRate, as the reference find_threshold does.
fun build(t: Trajectory, features: GestureFeatures = GestureFeatures.fit(t), threshold: Double? = null, candidates: Int = 24): VariableMarkovOracle {
val frames = features.frames(t)
if (threshold != null) return VariableMarkovOracle(t, features, frames, threshold)
return thresholdScan(t, features, frames, candidates).maxBy { it.second }.first
}fun thresholdScan(t: Trajectory, features: GestureFeatures = GestureFeatures.fit(t), candidates: Int = 24): List<Pair<VariableMarkovOracle, Double>>
Every candidate oracle with its information rate, for plotting the threshold choice.
fun thresholdScan(t: Trajectory, features: GestureFeatures = GestureFeatures.fit(t), candidates: Int = 24): List<Pair<VariableMarkovOracle, Double>> =
thresholdScan(t, features, features.frames(t), candidates)
private fun thresholdScan(t: Trajectory, features: GestureFeatures, frames: Array<DoubleArray>, candidates: Int): List<Pair<VariableMarkovOracle, Double>> {
// spread of typical distances, from a deterministic sample of pairs
val r = java.util.Random(0)
val ds = DoubleArray(min(4000, frames.size * 4)) { distance(frames[r.nextInt(frames.size)], frames[r.nextInt(frames.size)]) }.sorted()
val hi = ds[ds.size / 2].coerceAtLeast(1e-6)
val lo = hi / 200
return (0 until candidates).map { k ->
val th = lo * Math.pow(hi / lo, k / (candidates - 1.0))
val o = VariableMarkovOracle(t, features, frames, th)
o to o.totalInformationRate()
}
}
}
}data class Cue(val position: Int, val distance: Double, val context: Int)
A place in memory matching a cue: continue from sample position. Lower distance is better.
internal class IntArrayList
internal class IntArrayList {
private var a = IntArray(2)