function [fn_best, zeta_best, A_best, H_best, rep] = refine_modes_basin_mc_fast( ...
    f, Hc, fn0, zeta0, A0, varargin)
%REFINE_MODES_BASIN_MC_FAST  Faster basin-hop via decimation + quasi-random starts + early-stop.
% Serial-only implementation: no parfor; safe on base MATLAB.
%
% Inputs/Name-Values (superset of the original):
%   'FRF','Model','DCGain','Delay','Areal','MinSepPct','Trials',
%   'SigmaFnPct','SigmaZetaPct','WeightPeaks','PeakBand','PeakWeight',
%   'TemperatureScale','FitBaseline'
%
%   --- Speed knobs (all optional) ---
%   'Decimate'          (true)
%   'DecimTargetN'      (1200)
%   'FocusWinFactor'    (1.75)
%   'UseParallel'       (false)  % accepted but ignored (serial-only)
%   'Sobol'             (true)   % prefer Sobol; falls back to Halton; else rand
%   'EarlyStopPatience' (max(150, 0.25*Trials))
%   'RefineMaxIters'    (40)
%   'Seed'              ([])

% ---------------- options ----------------
p = inputParser;
p.addParameter('FRF','',@(s)ischar(s)||isstring(s));
p.addParameter('Model','auto',@(s)ischar(s)||isstring(s));
p.addParameter('DCGain',1+0i);
p.addParameter('Delay',0);
p.addParameter('Areal',false,@(b)islogical(b)&&isscalar(b));
p.addParameter('MinSepPct',2e-4,@(x)isnumeric(x)&&isscalar(x)&&x>=0);
p.addParameter('Trials',1200,@(x)isnumeric(x)&&isscalar(x)&&x>=0);
p.addParameter('SigmaFnPct',[0.12 0.01],@isnumeric);
p.addParameter('SigmaZetaPct',[0.25 0.03],@isnumeric);
p.addParameter('WeightPeaks',false,@(b)islogical(b)&&isscalar(b));
p.addParameter('PeakBand',[0.85 1.15],@isnumeric);
p.addParameter('PeakWeight',3.0,@(x)isnumeric(x)&&isscalar(x));
p.addParameter('TemperatureScale',0.01,@(x)isnumeric(x)&&isscalar(x));
p.addParameter('FitBaseline',false,@(b)islogical(b)&&isscalar(b));
% speed knobs
p.addParameter('Decimate',true,@(b)islogical(b)&&isscalar(b));
p.addParameter('DecimTargetN',1200,@(x)isnumeric(x)&&isscalar(x)&&x>=64);
p.addParameter('FocusWinFactor',1.75,@(x)isnumeric(x)&&isscalar(x)&&x>0);
p.addParameter('UseParallel',false,@(b)islogical(b)&&isscalar(b)); % accepted but ignored
p.addParameter('Sobol',true,@(b)islogical(b)&&isscalar(b));
p.addParameter('EarlyStopPatience',[],@(x)isnumeric(x)&&isscalar(x)&&x>=10);
p.addParameter('RefineMaxIters',40,@(x)isnumeric(x)&&isscalar(x)&&x>=5);
p.addParameter('Seed',[]);
p.parse(varargin{:});
o = p.Results;

% RNG
if ~isempty(o.Seed), rng(o.Seed,'twister'); end

f = double(f(:)); Hc = Hc(:); w = 2*pi*f;
Nm = numel(fn0); if isempty(A0), A0 = zeros(Nm,1); end; A0 = A0(:);

% ---- decimate grid for MC trials (keep focus around expected poles) ----
if o.Decimate && numel(f) > o.DecimTargetN
    Iwin = focus_windows(f, fn0, o.FocusWinFactor);
    Ifull = unique(Iwin(:));
    if numel(Ifull) > o.DecimTargetN
        stride = max(1, floor(numel(Ifull)/o.DecimTargetN));
        Ifull = Ifull(1:stride:end);
    end
    Idec = Ifull;
else
    Idec = (1:numel(f)).';
end
f_d  = f(Idec); w_d = w(Idec); H_d = Hc(Idec);

% ---- helpers ----
frfType = lower(char(o.FRF));
modelForm = lower(char(o.Model));
if strcmp(modelForm,'auto'), modelForm = iff(strcmp(frfType,''),'multiplicative','additive'); end
xfer_d  = @(A,fn,zeta) transfer_from_params_local_w(w_d, A, fn, zeta, frfType, modelForm, o.DCGain, o.Delay);
xfer    = @(A,fn,zeta) transfer_from_params_local_w(w,   A, fn, zeta, frfType, modelForm, o.DCGain, o.Delay);
crmse   = @(X,Y) sqrt(mean(abs(X-Y).^2));

% initial error on decimated grid
H0_d = xfer_d(A0, fn0, zeta0);
err_best_d = crmse(H_d, H0_d);
fn_best = fn0(:); zeta_best = zeta0(:); A_best = A0(:); H_best = xfer(A_best, fn_best, zeta_best);

% schedule
T = max(1,o.Trials);
sig_fn   = linspace(o.SigmaFnPct(1),   o.SigmaFnPct(end),   T);
sig_zeta = linspace(o.SigmaZetaPct(1), o.SigmaZetaPct(end), T);
useQMC   = o.Sobol && (exist('sobolset','file')==2 || exist('haltonset','file')==2);

if isempty(o.EarlyStopPatience), o.EarlyStopPatience = max(150, round(0.25*T)); end
no_improve = 0; best_so_far = err_best_d;

expDelay_d = (o.Delay~=0) * exp(1i*w_d*o.Delay) + (o.Delay==0)*1;

% ---- SERIAL trials with early stopping ----
for k = 1:T
    % quasi-Gaussian jitters with robust fallback
    z = qmc_gauss(2*Nm, k, useQMC);
    s_fn = sig_fn(k); s_zet = sig_zeta(k);
    d_fn  = s_fn  .* fn0(:).'            .* z(1:Nm);
    d_zet = s_zet .* max(1e-6,zeta0(:).').* z(Nm+1:end);

    fn_k   = fn0(:) + d_fn(:);
    zeta_k = zeta0(:) + d_zet(:);
    [fn_k, idx] = sort(fn_k); zeta_k = zeta_k(idx);
    fn_k(fn_k<=0) = min(fn0)*0.5; zeta_k = min(0.25, max(1e-5, zeta_k));

    % light separation enforcement
    for repSep = 1:2
        d = diff(fn_k);
        tooClose = d < (o.MinSepPct * 0.5) .* fn_k(1:end-1);
        if any(tooClose)
            I = find(tooClose);
            fn_k([I; I+1]) = fn_k([I; I+1]) + [-1; +1] .* 0.5 .* (o.MinSepPct .* fn_k([I; I+1]));
        else, break
        end
    end

    % LS residues on decimated grid
    Hcols = build_cols(w_d, fn_k, zeta_k, frfType);
    rhs   = build_rhs(H_d, o.DCGain, o.Delay, w_d, modelForm, expDelay_d);
    A_k   = Hcols \ rhs; if o.Areal, A_k = real(A_k); end

    % quick LM on decimated grid (if available)
    [fn_k, zeta_k, A_k] = quick_lm_decim(f_d, H_d, fn_k, zeta_k, A_k, o, frfType, modelForm);

    % error and accept
    H_est = synth_from_cols(Hcols, A_k, o.DCGain, modelForm, expDelay_d);
    err_k = sqrt(mean(abs(H_est - H_d).^2));

    if err_k < err_best_d*(1-1e-9)
        err_best_d = err_k; fn_best=fn_k; zeta_best=zeta_k; A_best=A_k;
        no_improve = 0; best_so_far = err_best_d;
    else
        no_improve = no_improve + 1;
        if no_improve >= o.EarlyStopPatience, break; end
    end
end
trials_run = k;

% ---- final full-grid LM polish for accuracy ----
[fn_best, zeta_best, A_best, H_best] = final_polish_fullgrid( ...
    f, Hc, fn_best, zeta_best, A_best, o);

rep = struct();
rep.decimated_N = numel(Idec);
rep.trials_run  = trials_run;
rep.rmse_decim_initial = best_so_far;
rep.rmse_decim_best    = err_best_d;
rep.rmse_after         = sqrt(mean(abs(Hc - H_best).^2));

% ------------- nested helpers -------------
function z = qmc_gauss(dim, k, useQMC)
    % Return a 1xDIM vector ~ N(0,1) via inverse-CDF of a quasi-random point.
    if useQMC
        try
            if exist('sobolset','file')==2
                ps = sobolset(dim); ps = scramble(ps,'MatousekAffineOwen');
                u = net(ps, k); u = u(end,:);
            elseif exist('haltonset','file')==2
                ps = haltonset(dim); ps = scramble(ps,'RR2');
                u = net(ps, k); u = u(end,:);
            else
                u = rand(1,dim);
            end
            u = min(max(u, realmin), 1-realmin); % guard
            z = sqrt(2) * erfinv(2*u - 1);
            return
        catch
            % fall through to rand
        end
    end
    z = randn(1,dim);
end

function Hcols = build_cols(wv, fn, zeta, frfType)
    Nd = numel(wv); Nm = numel(fn);
    Hcols = zeros(Nd, Nm);
    for m = 1:Nm
        wn = 2*pi*fn(m);
        den = (wn.^2 - wv.^2) + 1i*(2*zeta(m)*wn).*wv;
        switch frfType
            case 'disp/force', Hcols(:,m) = 1 ./ den;
            case 'vel/force',  Hcols(:,m) = (1i*wv) ./ den;
            case 'acc/force',  Hcols(:,m) = (-wv.^2) ./ den;
            otherwise,         Hcols(:,m) = (wv.^2) ./ den;
        end
    end
end

function rhs = build_rhs(Hy, G, tau, wv, modelForm, expDelay_v)
    if tau~=0, Hy = Hy .* expDelay_v; end
    if strcmp(modelForm,'multiplicative') && G~=0, rhs = Hy./G - 1; else, rhs = Hy - G; end
end

function Hsyn = synth_from_cols(Hcols, A, G, modelForm, expDelay_v)
    body = Hcols*A;
    if strcmp(modelForm,'multiplicative'), Hsyn = G * (1 + body); else, Hsyn = G + body; end
    Hsyn = Hsyn ./ expDelay_v; % undo delay factor used in rhs
end

function [fn_q, zeta_q, A_q] = quick_lm_decim(f_d, H_d, fn_i, zeta_i, A_i, o, frfType, modelForm)
    fn_q = fn_i; zeta_q = zeta_i; A_q = A_i;
    if exist('refine_modes_levenberg','file')==2
        [fn_q, zeta_q, A_q] = refine_modes_levenberg( ...
            f_d, H_d, fn_i, zeta_i, A_i, ...
            'FRF', frfType, 'Model', modelForm, 'DCGain', o.DCGain, ...
            'Areal', false, 'Delay', o.Delay, 'EstimateDelay', false, ...
            'FitBaseline', false, 'MaxIters', o.RefineMaxIters, 'Lambda0', 1e-2, ...
            'MinSepPct', o.MinSepPct, 'Verbose', false); %#ok<ASGLU>
    end
end

function Hhat = transfer_from_params_local_w(wv, A, fn, zeta, frfType, modelForm, G, tau)
    Nd = numel(wv); Nm = numel(fn);
    Hcols = zeros(Nd, Nm);
    for m = 1:Nm
        wn = 2*pi*fn(m);
        den = (wn.^2 - wv.^2) + 1i*(2*zeta(m)*wn).*wv;
        switch frfType
            case 'disp/force', Hcols(:,m) = 1 ./ den;
            case 'vel/force',  Hcols(:,m) = (1i*wv) ./ den;
            case 'acc/force',  Hcols(:,m) = (-wv.^2) ./ den;
            otherwise,         Hcols(:,m) = (wv.^2) ./ den;
        end
    end
    body = Hcols*A;
    if strcmp(modelForm,'multiplicative'), Hhat = G * (1 + body); else, Hhat = G + body; end
    if tau~=0, Hhat = exp(-1i*wv*tau).*Hhat; end
end

function [fn_f, zeta_f, A_f, H_f] = final_polish_fullgrid(f, Hc, fn_b, zeta_b, A_b, o)
    if exist('refine_modes_levenberg','file')==2
        [fn_f, zeta_f, A_f, H_f] = refine_modes_levenberg( ...
            f, Hc, fn_b, zeta_b, A_b, ...
            'FRF', o.FRF, 'Model', o.Model, 'DCGain', o.DCGain, ...
            'Areal', false, 'Delay', o.Delay, 'EstimateDelay', true, ...
            'FitBaseline', false, 'MaxIters', 120, 'Lambda0', 1e-2, ...
            'MinSepPct', o.MinSepPct, 'Verbose', false); %#ok<ASGLU>
    else
        A_f = A_b; fn_f = fn_b; zeta_f = zeta_b;
        H_f = transfer_from_params_local_w(2*pi*f, A_f, fn_f, zeta_f, lower(char(o.FRF)), lower(char(o.Model)), o.DCGain, o.Delay);
    end
end

function I = focus_windows(fv, fnc, winFact)
    I = [];
    for m = 1:numel(fnc)
        f0 = fnc(m);
        fLo = max(fv(1),  f0*(1 - 1/winFact));
        fHi = min(fv(end),f0*(1 + 1/winFact));
        I = [I; find(fv>=fLo & fv<=fHi)]; %#ok<AGROW>
    end
    if isempty(I), I = (1:numel(fv)).'; end
    I = unique(I(:));
end

function v = iff(c,a,b), if c, v=a; else, v=b; end, end
end
