function [fcand, zcand] = preidentify_modes_complex_local(f, Hc, kmax, varargin)
% PREIDENTIFY_MODES_COMPLEX_LOCAL
% Three cues: |Im(H)|, phase curvature, |H-1|.
% 1) Find raw peaks
% 2) Consolidate near-duplicates into frequency clusters
% 3) FORCE exactly kmax seeds via separation-aware selection
% 4) Rough zeta from |H-1| half-height bandwidth

p = inputParser;
p.addParameter('MinSepPct',       3e-4);   % used inside peak finders
p.addParameter('ConsolidatePct',  6e-3);   % ~0.6% cluster width
p.addParameter('Verbose',         false);
p.parse(varargin{:});
opt = p.Results;

% ---- sanitize inputs ----
f  = double(f(:));
Hc = Hc(:);
N  = numel(f);
df = max(eps, median(diff(f)));
dmin = max(3, round(opt.MinSepPct * N));  % samples

% ---------- CUE 1: |Im(H)| ----------
yI = abs(imag(Hc));
yI = smooth1d_local(yI, max(5, round(0.002*N)));
[pI, lI] = findpeaks_local(yI, ...
    'MinPeakProminence', max(1e-3*max(yI), eps), ...
    'MinPeakDistance',   dmin);

% ---------- CUE 2: phase curvature ----------
ph = unwrap(angle(Hc));
d2 = abs([0; diff(ph,2); 0]);
d2 = smooth1d_local(d2, max(5, round(0.0015*N)));
[pP, lP] = findpeaks_local(d2, ...
    'MinPeakProminence', max(1e-3*max(d2), eps), ...
    'MinPeakDistance',   dmin);

% ---------- CUE 3: |H-1| ----------
R  = abs(Hc - 1);
Rm = max(R); if Rm==0 || ~isfinite(Rm), Rm=1; end
R  = R ./ Rm;
R  = smooth1d_local(R, max(3, round(0.0012*N)));
dR = diff(R); meddR = 0; if ~isempty(dR), meddR = median(abs(dR)); end
promR = max(0.05*max(R), 5*meddR + 1e-6);
[pR, lR] = findpeaks_local(R, ...
    'MinPeakProminence', promR, ...
    'MinPeakDistance',   max(dmin, round(0.002*N)));

% ---------- union + fused score ----------
Lall = unique([lI(:); lP(:); lR(:)]);
if isempty(Lall), fcand=[]; zcand=[]; return; end
fi = f(Lall);

sI = zeros(size(Lall)); sP = sI; sR = sI;
if ~isempty(lI), sI = interp1(lI, pI, Lall, 'nearest', 0); end
if ~isempty(lP), sP = interp1(lP, pP, Lall, 'nearest', 0); end
if ~isempty(lR), sR = interp1(lR, pR, Lall, 'nearest', 0); end
if any(sI), sI = sI./max(sI); end
if any(sP), sP = sP./max(sP); end
if any(sR), sR = sR./max(sR); end
wI = 1.00; wP = 0.70; wR = 1.15;
score = wI*sI + wP*sP + wR*sR;

% ---------- consolidate near-duplicates into clusters ----------
[fs, ord] = sort(fi); ss = score(ord);
clusters = {};
s = 1;
for i = 2:numel(fs)
    if (fs(i)-fs(i-1)) <= opt.ConsolidatePct*max(fs(i),fs(i-1))
        % same cluster
    else
        clusters{end+1} = s:(i-1); %#ok<AGROW>
        s = i;
    end
end
clusters{end+1} = s:numel(fs); %#ok<AGROW>

% representative of each cluster = highest fused score in the cluster
nc = numel(clusters);
rep_f = zeros(nc,1);
rep_s = zeros(nc,1);
for c = 1:nc
    idx  = clusters{c};
    [~, jloc] = max(ss(idx));
    rep_f(c) = fs(idx(jloc));
    rep_s(c) = ss(idx(jloc));
end

% ---------- FORCE exactly kmax seeds ----------
if nc <= kmax
    sel_f = rep_f;
    sel_s = rep_s;
else
    % Seed with best-scoring rep, then farthest-point sampling to maximize separation
    [~, j0] = max(rep_s);
    sel = j0;
    while numel(sel) < kmax
        % distance (in Hz) to nearest already-selected rep
        d = inf(nc,1);
        for j = 1:nc
            d(j) = min(abs(rep_f(j) - rep_f(sel)));
        end
        d(sel) = -inf;           % exclude already selected
        % small score bonus to break ties toward stronger peaks
        [~, jnext] = max(d + 0.01*rep_s);
        sel(end+1) = jnext; %#ok<AGROW>
    end
    sel = unique(sel, 'stable');
    sel_f = rep_f(sel);
    sel_s = rep_s(sel); %#ok<NASGU>
end

fc = sel_f(:);
if isempty(fc), fcand=[]; zcand=[]; return; end

% ---------- estimate ζ from |H-1| half-height (robust, no interp1) ----------
zc = 0.02*ones(size(fc));
for k = 1:numel(fc)
    [~, i0] = min(abs(f - fc(k)));
    pk = R(i0); if ~(isfinite(pk) && pk>0), continue; end
    hh = 0.5 * pk;

    % search left/right for the half-height crossings
    iL=i0; while iL>1        && R(iL) >= hh, iL=iL-1; end
    iR=i0; while iR<numel(R) && R(iR) >= hh, iR=iR+1; end

    if iL>1 && iR<numel(R)
        % linear two-point solve (avoid interp1 monotonicity constraints)
        fL = f(iL)   + (hh - R(iL))   * (f(iL+1)-f(iL))   / (R(iL+1)-R(iL)   + eps);
        fR = f(iR-1) + (hh - R(iR-1)) * (f(iR)-f(iR-1))   / (R(iR)-R(iR-1)   + eps);
        bw = max(df, fR - fL);
        zc(k) = max(5e-3, min(0.12, bw/(2*max(fc(k), eps))));
    end
end

[fcand, o2] = sort(fc(:));
zcand = zc(o2);

if opt.Verbose
    fprintf('[pre-ID] seeds (forced %d):', numel(fcand));
    fprintf(' %.1f', fcand); fprintf('\n');
end
end

% --- local helpers ---
function y = smooth1d_local(y, win)
win = max(3, round(win)); win = min(win, numel(y));
k = ones(win,1)/max(1,win);
y = conv(double(y(:)), k, 'same');
end

function [pks, locs] = findpeaks_local(y, varargin)
if exist('findpeaks','file')==2
    [pks,locs] = findpeaks(y, varargin{:});
    return
end
y = y(:); n = numel(y);
locs = find([false; (y(2:n-1) > y(1:n-2)) & (y(2:n-1) >= y(3:n)); false]);
pks  = y(locs);
MPP=[]; MPD=[];
for k=1:2:numel(varargin)
    if strcmpi(varargin{k},'MinPeakProminence'), MPP=varargin{k+1}; end
    if strcmpi(varargin{k},'MinPeakDistance'),   MPD=varargin{k+1}; end
end
if ~isempty(MPP)
    keep = pks >= MPP;
    pks = pks(keep); locs = locs(keep);
end
if ~isempty(MPD) && numel(locs)>1
    keep = true(size(locs));
    last = -Inf;
    for i=1:numel(locs)
        if locs(i) - last < MPD, keep(i)=false; else, last = locs(i); end
    end
    pks = pks(keep); locs = locs(keep);
end
end
