function [times, labels] = overlap_pursuit(x, fs, training)
% Silence segmentation plus per-block matching pursuit for overlaps.
%
% The recordings place each event, an exact snippet scaled by one gain,
% on a digitally silent background, so activity splits into blocks: a
% block exactly one snippet long holds a single event, and a longer
% block holds several overlapping events whose first onset is the block
% start and whose last onset is the block end minus one snippet length.
% Singles are classified by the spectral kNN of silence_gap_knn. For a
% multi-event block the two boundary onsets are taken as anchors and
% fitted by least squares against every training snippet (the best
% instance and gain per onset, refined by a few coordinate-descent
% sweeps that add each event back and refit it against the residual).
% Two pursuit stages then add interior events. First, coverage: every
% sample in a block must lie inside some event's snippet span, so while
% uncovered samples remain, the residual is scanned by FFT
% cross-correlation against the whole bank and the best-reducing onset
% whose span touches uncovered samples is added; this is a proof an
% interior event exists, not a guess. Second, duplicates: remaining
% residual energy is scanned the same way, and a candidate is accepted
% while its gain and energy reduction are plausible, unless an event of
% the same type already sits within a few hundred samples, since
% instance-to-instance mismatch of an already-fitted event is the usual
% source of such candidates and near-coincident same-type pairs are
% rare in comparison.
%
% Limitations: the block bookkeeping relies on the silent background,
% and nearly-coincident events of the same type are reported once. The
% snippet bank explains an event only up to instance variability
% (worst for the least stereotyped type), which is what bounds how far
% the residual can be trusted; with added background noise the
% thresholds below would need retuning.
nTypes = numel(training.snips);
L = size(training.snips{1}, 1);
pre = training.pre;
n = numel(x);

% tuned constants
gmin = 0.35;        % smallest plausible fitted gain for a pursuit event
redmin = 0.5;       % smallest energy reduction for a duplicate event
prox = round(0.0156 * fs);   % same-type exclusion radius, ~250 samples
maxExtra = 6;       % coverage-pursuit iterations per block
maxDup = 4;         % duplicate-pursuit iterations per block

% ---- training bank ------------------------------------------------------
nTrain = 0;
for k = 1:nTypes
  nTrain = nTrain + size(training.snips{k}, 2);
end
B = zeros(L, nTrain);
bt = zeros(nTrain, 1);
c = 0;
for k = 1:nTypes
  S = training.snips{k};
  for j = 1:size(S, 2)
    c = c + 1;
    B(:, c) = S(:, j);
    bt(c) = k;
  end
end
n2 = zeros(nTrain, 1);
for j = 1:nTrain
  n2(j) = sum(B(:, j) .* B(:, j));
end

% ---- kNN training features for single-event blocks ----------------------
F = zeros(nTrain, 8);
for j = 1:nTrain
  F(j, :) = spectral_features(B(pre+1:end, j), fs);
end
mu = mean(F, 1);
sd = std(F, 0, 1) + 1e-9;
Fz = (F - repmat(mu, nTrain, 1)) ./ repmat(sd, nTrain, 1);

% ---- silence segmentation ------------------------------------------------
eps0 = 1e-4 * max(abs(x));
gapClose = round(0.002 * fs);
minDur = round(0.010 * fs);
bs = zeros(0, 1);
be = zeros(0, 1);
i = 1;
while i <= n
  if abs(x(i)) > eps0
    s = i;
    e = i;
    j = i + 1;
    quiet = 0;
    while j <= n && quiet <= gapClose
      if abs(x(j)) > eps0
        e = j;
        quiet = 0;
      else
        quiet = quiet + 1;
      end
      j = j + 1;
    end
    if e - s + 1 >= minDur
      bs(end+1, 1) = s;
      be(end+1, 1) = e;
    end
    i = e + gapClose + 1;
  else
    i = i + 1;
  end
end

% ---- per-block detection and classification ------------------------------
times = zeros(0, 1);
labels = zeros(0, 1);
for bi = 1:numel(bs)
  s = bs(bi);
  D = be(bi) - s + 1;
  if D <= L + 2
    % single event: onset at the block start, spectral kNN label
    seg = x(s+pre : min(n, s+L-1));
    f = (spectral_features(seg, fs) - mu) ./ sd;
    d = sqrt(sum((Fz - repmat(f, nTrain, 1)).^2, 2));
    [ds, di] = sort(d);
    kNN = min(5, nTrain);
    votes = zeros(nTypes, 1);
    for v = 1:kNN
      votes(bt(di(v))) = votes(bt(di(v))) + 1 / (ds(v) + 1e-6);
    end
    [~, lab] = max(votes);
    times(end+1, 1) = (s - 1 + pre) / fs;
    labels(end+1, 1) = lab;
  else
    blk = x(s : s+D-1);
    [evOn, evType] = decompose(blk, B, bt, n2, L, gmin, redmin, prox, maxExtra, maxDup);
    for q = 1:numel(evOn)
      times(end+1, 1) = (s - 1 + evOn(q) - 1 + pre) / fs;
      labels(end+1, 1) = evType(q);
    end
  end
end
end

function [evOn, evType] = decompose(blk, B, bt, n2, L, gmin, redmin, prox, maxExtra, maxDup)
% Anchored matching pursuit over one multi-event block.
D = numel(blk);
nl = D - L + 1;
nTrain = numel(bt);

% conjugate template spectra for the FFT correlation scans
nfft = 2^nextpow2(D);
BF = zeros(nfft, nTrain) + 0i;
for j = 1:nTrain
  BF(:, j) = conj(fft(B(:, j), nfft));
end

% events as parallel arrays, seeded with the two exact boundary anchors
evOn = [1; nl];
evType = zeros(2, 1);
evIdx = zeros(2, 1);
evGain = zeros(2, 1);
residual = blk;

[residual, evType, evIdx, evGain] = sweep_fit(residual, evOn, evType, evIdx, evGain, B, bt, n2, L, 2);

% coverage pursuit: every block sample must lie in some event span
for it = 1:maxExtra
  cov = zeros(D, 1);
  for q = 1:numel(evOn)
    cov(evOn(q) : evOn(q)+L-1) = 1;
  end
  nUnc = D - sum(cov);
  if nUnc == 0
    break;
  end
  % lags whose span [o, o+L-1] contains an uncovered sample
  m = conv(1 - cov, ones(L, 1));
  mask = zeros(nl, 1);
  for o = 1:nl
    mask(o) = m(o + L - 1) > 0;
  end
  [o, kk, idx, g, red] = scan_best(residual, BF, bt, n2, L, mask);
  if idx == 0
    break;
  end
  residual(o:o+L-1) = residual(o:o+L-1) - g * B(:, idx);
  evOn(end+1, 1) = o;
  evType(end+1, 1) = kk;
  evIdx(end+1, 1) = idx;
  evGain(end+1, 1) = g;
  [residual, evType, evIdx, evGain] = sweep_fit(residual, evOn, evType, evIdx, evGain, B, bt, n2, L, 2);
end

% duplicate pursuit: events fully covered by other spans
mask = zeros(nl, 1) + 1;
for it = 1:maxDup
  [o, kk, idx, g, red] = scan_best(residual, BF, bt, n2, L, mask);
  if idx == 0 || red < redmin || g < gmin
    break;
  end
  clash = 0;
  for q = 1:numel(evOn)
    if evOn(q) == o || (abs(evOn(q) - o) < prox && evType(q) == kk)
      clash = 1;
    end
  end
  if clash
    break;
  end
  residual(o:o+L-1) = residual(o:o+L-1) - g * B(:, idx);
  evOn(end+1, 1) = o;
  evType(end+1, 1) = kk;
  evIdx(end+1, 1) = idx;
  evGain(end+1, 1) = g;
  [residual, evType, evIdx, evGain] = sweep_fit(residual, evOn, evType, evIdx, evGain, B, bt, n2, L, 2);
end
end

function [residual, evType, evIdx, evGain] = sweep_fit(residual, evOn, evType, evIdx, evGain, B, bt, n2, L, nsw)
% Coordinate descent: add each event back and refit its type/instance/gain.
[~, ord] = sort(evOn);
for it = 1:nsw
  for z = 1:numel(ord)
    q = ord(z);
    o = evOn(q);
    if evIdx(q) > 0
      residual(o:o+L-1) = residual(o:o+L-1) + evGain(q) * B(:, evIdx(q));
    end
    seg = residual(o:o+L-1);
    bestRed = 0;
    bestJ = 0;
    bestG = 0;
    for j = 1:numel(bt)
      cj = sum(seg .* B(:, j));
      if cj > 0
        red = cj * cj / n2(j);
        if red > bestRed
          bestRed = red;
          bestJ = j;
          bestG = cj / n2(j);
        end
      end
    end
    if bestJ > 0
      residual(o:o+L-1) = residual(o:o+L-1) - bestG * B(:, bestJ);
      evType(q) = bt(bestJ);
      evIdx(q) = bestJ;
      evGain(q) = bestG;
    end
  end
end
end

function [o, kk, idx, g, red] = scan_best(residual, BF, bt, n2, L, mask)
% Best (onset, template) over masked lags by FFT cross-correlation;
% reduction is the residual energy removed by the least-squares gain.
nfft = size(BF, 1);
nl = numel(mask);
R = fft(residual, nfft);
o = 0; kk = 0; idx = 0; g = 0; red = 0;
for j = 1:numel(bt)
  cf = real(ifft(R .* BF(:, j)));
  cl = cf(1:nl);
  rv = (cl .* cl) / n2(j);
  rv = rv .* (cl > 0) .* mask;
  [m, mi] = max(rv);
  if m > red
    red = m;
    o = mi;
    kk = bt(j);
    idx = j;
    g = cl(mi) / n2(j);
  end
end
end

function f = spectral_features(seg, fs)
% Seven log band-energy ratios (mean removed, level cancels) and the
% log spectral centroid; the extra low bands split the two bass types.
nfft = 2048;
m = min(numel(seg), nfft);
s2 = seg(1:m);
win = 0.5 - 0.5 * cos(2 * pi * (0:m-1)' / max(m - 1, 1));
S = fft(s2 .* win, nfft);
p = abs(S(1:nfft/2+1)).^2;
hz = fs / nfft;
edges = [0 100 200 400 800 1600 3200 6400];
b = zeros(1, 7);
for q = 1:7
  k0 = max(2, floor(edges(q) / hz) + 1);
  k1 = min(nfft/2 + 1, floor(edges(q+1) / hz));
  b(q) = log10(sum(p(k0:k1)) + 1e-12);
end
b = b - mean(b);
kv = (0:numel(p)-1)' * hz;
cen = sum(p .* kv) / (sum(p) + 1e-12);
f = [b, log10(cen + 1)];
end