| 1 | import { CATEGORIES, clamp, errorOnsetOf, quantile, type BinLevel, type Metric, type SignalPyramid, type WhaleEvent } from './model.js'; |
| 2 | |
| 3 | function emptyLevel(length: number, binMs: number): BinLevel { |
| 4 | const n = CATEGORIES.length * length; |
| 5 | return { binMs, length, onsets: new Float64Array(n), activeMs: new Float64Array(n), |
| 6 | outputTokens: new Float64Array(n), cost: new Float64Array(n), errors: new Float64Array(n), peak: new Float64Array(n) }; |
| 7 | } |
| 8 | const FIELDS = ['onsets', 'activeMs', 'outputTokens', 'cost', 'errors'] as const; |
| 9 | /** O(events + channels × bins), including long intervals. No span-length inner loop. */ |
| 10 | export function buildPyramid(events: readonly WhaleEvent[], requestedDuration?: number, maxBins = 16_384): SignalPyramid { |
| 11 | if (!Number.isInteger(maxBins) || maxBins < 16 || maxBins > 1_048_576) throw new Error('maxBins must be an integer in [16, 1048576].'); |
| 12 | let duration = requestedDuration ?? 1; |
| 13 | for (const e of events) duration = Math.max(duration, e.endTime, e.startTime, e.status === 'error' ? errorOnsetOf(e) : 0); |
| 14 | if (!Number.isFinite(duration) || duration < 0) throw new Error('Signal duration must be finite and nonnegative.'); |
| 15 | duration = Math.max(1, duration); |
| 16 | const binMs = 2 ** Math.ceil(Math.log2(Math.max(1, duration / (maxBins - 1)))); |
| 17 | const length = Math.floor(duration / binMs) + 1, fine = emptyLevel(length, binMs); |
| 18 | const stride = length + 1; |
| 19 | const activeDiff = new Float64Array(CATEGORIES.length * stride), tokenDiff = new Float64Array(CATEGORIES.length * stride); |
| 20 | for (const e of events) { |
| 21 | if (!Number.isFinite(e.startTime) || !Number.isFinite(e.endTime) || e.startTime < 0 || e.endTime < e.startTime) throw new Error(`Invalid interval for ${e.id}.`); |
| 22 | const channel = CATEGORIES.indexOf(e.category); |
| 23 | if (channel < 0) throw new Error(`Unknown category for ${e.id}.`); |
| 24 | const a = Math.floor(e.startTime / binMs), b = Math.floor(e.endTime / binMs); |
| 25 | const at = channel * length, diff = channel * stride; |
| 26 | fine.onsets[at + a]++; |
| 27 | fine.cost[at + a] += e.cost ?? 0; |
| 28 | if (e.status === 'error') fine.errors[at + Math.floor(errorOnsetOf(e) / binMs)]++; |
| 29 | const d = e.endTime - e.startTime, tokens = e.outputTokens ?? 0; |
| 30 | if (d === 0) { fine.outputTokens[at + a] += tokens; continue; } |
| 31 | if (a === b) { fine.activeMs[at + a] += d; fine.outputTokens[at + a] += tokens; } |
| 32 | else { |
| 33 | const left = (a + 1) * binMs - e.startTime, right = e.endTime - b * binMs, tokenRate = tokens / d; |
| 34 | fine.activeMs[at + a] += left; |
| 35 | fine.activeMs[at + b] += right; |
| 36 | fine.outputTokens[at + a] += left * tokenRate; |
| 37 | fine.outputTokens[at + b] += right * tokenRate; |
| 38 | if (b > a + 1) { |
| 39 | activeDiff[diff + a + 1] += binMs; activeDiff[diff + b] -= binMs; |
| 40 | tokenDiff[diff + a + 1] += binMs * tokenRate; tokenDiff[diff + b] -= binMs * tokenRate; |
| 41 | } |
| 42 | } |
| 43 | } |
| 44 | for (let c = 0; c < CATEGORIES.length; c++) { |
| 45 | let active = 0, tokens = 0; |
| 46 | for (let i = 0; i < length; i++) { |
| 47 | active += activeDiff[c * stride + i]; tokens += tokenDiff[c * stride + i]; |
| 48 | const at = c * length + i; |
| 49 | fine.activeMs[at] = Math.max(0, fine.activeMs[at] + active); |
| 50 | fine.outputTokens[at] = Math.max(0, fine.outputTokens[at] + tokens); |
| 51 | fine.peak[at] = fine.activeMs[at] / binMs; |
| 52 | } |
| 53 | } |
| 54 | const levels = [fine]; |
| 55 | while (levels.at(-1)!.length > 1) { |
| 56 | const child = levels.at(-1)!, parent = emptyLevel(Math.ceil(child.length / 2), child.binMs * 2); |
| 57 | for (let c = 0; c < CATEGORIES.length; c++) for (let i = 0; i < parent.length; i++) { |
| 58 | const a = c * child.length + i * 2, b = a + 1, dst = c * parent.length + i, hasB = i * 2 + 1 < child.length; |
| 59 | for (const key of FIELDS) parent[key][dst] = child[key][a] + (hasB ? child[key][b] : 0); |
| 60 | parent.peak[dst] = Math.max(child.peak[a], hasB ? child.peak[b] : 0); |
| 61 | } |
| 62 | levels.push(parent); |
| 63 | } |
| 64 | const p: SignalPyramid = { duration, channels: CATEGORIES, levels, calibration: { activity: 1, onsets: 1, tokens: 1, cost: 1 } }; |
| 65 | const reference = chooseLevel(p, duration / 1200); |
| 66 | for (const metric of ['activity', 'onsets', 'tokens', 'cost'] as Metric[]) { |
| 67 | const positives: number[] = []; |
| 68 | for (let c = 0; c < CATEGORIES.length; c++) for (let i = 0; i < reference.length; i++) { |
| 69 | const v = binValue(reference, c, i, metric); if (v > 0) positives.push(v); |
| 70 | } |
| 71 | p.calibration[metric] = Math.max(1e-9, quantile(positives, .95)); |
| 72 | } |
| 73 | return p; |
| 74 | } |
| 75 | /** Choose the coarsest stored level no wider than one requested pixel interval. */ |
| 76 | export function chooseLevel(p: SignalPyramid, targetBinMs: number): BinLevel { |
| 77 | let result = p.levels[0]; |
| 78 | for (const level of p.levels) { if (level.binMs > targetBinMs) break; result = level; } |
| 79 | return result; |
| 80 | } |
| 81 | export function binValue(level: BinLevel, channel: number, bin: number, metric: Metric): number { |
| 82 | if (bin < 0 || bin >= level.length) return 0; |
| 83 | const at = channel * level.length + bin; |
| 84 | if (metric === 'activity') return level.activeMs[at] / level.binMs; |
| 85 | if (metric === 'onsets') return level.onsets[at] * 1000 / level.binMs; |
| 86 | if (metric === 'tokens') return level.outputTokens[at] * 1000 / level.binMs; |
| 87 | return level.cost[at] * 1000 / level.binMs; |
| 88 | } |
| 89 | export function intensity(value: number, reference: number): number { |
| 90 | return clamp(Math.log1p(value / Math.max(1e-9, reference) * 8) / Math.log(9), 0, 1); |
| 91 | } |
| 92 | /** Conserved totals: useful for tests and alternate native backends. */ |
| 93 | export function totals(level: BinLevel): Record<string, number> { |
| 94 | return Object.fromEntries(FIELDS.map(k => [k, level[k].reduce((s, x) => s + x, 0)])); |
| 95 | } |
| 96 | |
| 97 | /** Sorted-start, max-end segment tree. Long root spans do not force a reverse |
| 98 | * scan through every earlier event. Results are chronological and honor limits. */ |
| 99 | export class IntervalIndex { |
| 100 | readonly events: WhaleEvent[]; |
| 101 | private readonly maxEnd: Float64Array; |
| 102 | private readonly leafCount: number; |
| 103 | constructor(events: readonly WhaleEvent[]) { |
| 104 | this.events = [...events].sort((a, b) => a.startTime - b.startTime || a.id.localeCompare(b.id)); |
| 105 | this.leafCount = 2 ** Math.ceil(Math.log2(Math.max(1, this.events.length))); |
| 106 | this.maxEnd = new Float64Array(this.leafCount * 2).fill(-Infinity); |
| 107 | for (let i = 0; i < this.events.length; i++) this.maxEnd[this.leafCount + i] = this.events[i].endTime; |
| 108 | for (let i = this.leafCount - 1; i; i--) this.maxEnd[i] = Math.max(this.maxEnd[i * 2], this.maxEnd[i * 2 + 1]); |
| 109 | } |
| 110 | query(start: number, end: number, limit = Infinity): WhaleEvent[] { |
| 111 | if (!Number.isFinite(start) || !Number.isFinite(end) || end < start || limit <= 0) return []; |
| 112 | let lo = 0, hi = this.events.length; |
| 113 | while (lo < hi) { const mid = (lo + hi) >>> 1; if (this.events[mid].startTime <= end) lo = mid + 1; else hi = mid; } |
| 114 | const bound = lo, found: WhaleEvent[] = []; |
| 115 | const visit = (node: number, left: number, right: number): void => { |
| 116 | if (left >= bound || this.maxEnd[node] < start || found.length >= limit) return; |
| 117 | if (right - left === 1) { |
| 118 | const e = this.events[left]; |
| 119 | if (e && (e.endTime > start || e.startTime === e.endTime && e.startTime >= start)) found.push(e); |
| 120 | return; |
| 121 | } |
| 122 | const mid = (left + right) >>> 1; |
| 123 | visit(node * 2, left, mid); visit(node * 2 + 1, mid, right); |
| 124 | }; |
| 125 | visit(1, 0, this.leafCount); |
| 126 | return found; |
| 127 | } |
| 128 | } |
| 129 | /** Sample an onset train into equal-width bins; nothing is inferred between impulses. */ |
| 130 | export function onsetSeries(events: readonly WhaleEvent[], start: number, end: number, n = 256, category?: string): Float64Array { |
| 131 | const out = new Float64Array(n), span = Math.max(1e-6, end - start); |
| 132 | for (const e of events) if ((!category || e.category === category) && e.startTime >= start && e.startTime < end) { |
| 133 | const at = Math.floor((e.startTime - start) / span * n); if (at >= 0 && at < n) out[at]++; |
| 134 | } |
| 135 | return out; |
| 136 | } |
| 137 | /** Centered, variance-normalized autocorrelation; lag-zero is one unless constant. */ |
| 138 | export function autocorrelation(input: ArrayLike<number>, maxLag = 64): number[] { |
| 139 | const n = input.length; |
| 140 | if (!n) return []; |
| 141 | let mean = 0; for (let i = 0; i < n; i++) mean += input[i]; mean /= n; |
| 142 | const centered = Array.from(input, x => x - mean), energy = centered.reduce((s, v) => s + v * v, 0); |
| 143 | const result: number[] = []; |
| 144 | for (let lag = 0; lag <= Math.min(maxLag, n - 1); lag++) { |
| 145 | let sum = 0; for (let i = 0; i < n - lag; i++) sum += centered[i] * centered[i + lag]; |
| 146 | result.push(energy > 1e-12 ? sum / energy : 0); |
| 147 | } |
| 148 | return result; |
| 149 | } |
| 150 | export interface Spectrum { frequencies: number[]; power: number[]; entropy: number; peakHz: number; sampleHz: number; resolutionHz: number } |
| 151 | /** Hann-window periodogram of an actual uniformly sampled onset signal, DC removed. */ |
| 152 | export function periodogram(input: ArrayLike<number>, sampleHz: number): Spectrum { |
| 153 | const n = input.length, frequencies: number[] = [], power: number[] = []; |
| 154 | if (n < 4 || sampleHz <= 0) return { frequencies, power, entropy: 0, peakHz: 0, sampleHz, resolutionHz: 0 }; |
| 155 | let mean = 0; for (let i = 0; i < n; i++) mean += input[i]; mean /= n; |
| 156 | const windowed = Array.from(input, (x, i) => (x - mean) * (.5 - .5 * Math.cos(2 * Math.PI * i / (n - 1)))); |
| 157 | const windowEnergy=Array.from({length:n},(_,i)=>(.5-.5*Math.cos(2*Math.PI*i/(n-1)))**2).reduce((a,b)=>a+b,0); |
| 158 | for (let k = 1; k <= Math.floor(n / 2); k++) { |
| 159 | let re = 0, im = 0; |
| 160 | for (let t = 0; t < n; t++) { const angle = 2 * Math.PI * k * t / n; re += windowed[t] * Math.cos(angle); im -= windowed[t] * Math.sin(angle); } |
| 161 | frequencies.push(k * sampleHz / n); power.push((re * re + im * im) / (windowEnergy * sampleHz) * (n%2===0&&k===n/2?1:2)); |
| 162 | } |
| 163 | const sum = power.reduce((s, x) => s + x, 0); |
| 164 | let entropy = 0; for (const x of power) if (x > 0 && sum > 0) { const p = x / sum; entropy -= p * Math.log2(p); } |
| 165 | entropy = power.length > 1 ? entropy / Math.log2(power.length) : 0; |
| 166 | const max = Math.max(...power); |
| 167 | return { frequencies, power, entropy, peakHz: max > 1e-12 ? frequencies[power.indexOf(max)] : 0, sampleHz, resolutionHz: sampleHz / n }; |
| 168 | } |
| 169 | export function unionDuration(events: readonly WhaleEvent[], start = 0, end = Infinity): number { |
| 170 | const intervals = events.filter(e => e.endTime > e.startTime && e.endTime > start && e.startTime < end) |
| 171 | .map(e => [Math.max(start, e.startTime), Math.min(end, e.endTime)]).sort((a, b) => a[0] - b[0]); |
| 172 | let total = 0, left = 0, right = 0; |
| 173 | for (const [a, b] of intervals) { |
| 174 | if (a > right) { total += right - left; left = a; right = b; } else right = Math.max(right, b); |
| 175 | } |
| 176 | return total + right - left; |
| 177 | } |
| 178 |