function [fn, zeta, A, Hhat, report] = refine_modes_basin_mc( ...
        f, Hc, fn0, zeta0, A0, varargin)
% REFINE_MODES_BASIN_MC  Monte-Carlo "basin hop" around an initial pole set.
%   [fn, zeta, A, Hhat, report] = refine_modes_basin_mc(f, Hc, fn0, zeta0, A0, ...)
%
% Required:
%   f       : frequency vector (Hz), column or row
%   Hc      : complex FRF (same length as f), column or row
%   fn0     : initial natural frequencies (Hz), column
%   zeta0   : initial damping ratios, column
%   A0      : initial residues (complex unless 'Areal' = true), column
%
% Name-Value options (key ones):
%   'FRF'              '', 'acc/force', 'vel/force', 'disp/force'   (default '')
%   'Model'            'auto'|'additive'|'multiplicative'           (default 'auto')
%   'DCGain'           complex scalar (default 1 for '', 0 otherwise)
%   'Delay'            time delay [s] (default 0)
%   'Areal'            logical, force real residues                  (default false)
%   'MinSepPct'        minimum relative frequency separation         (default 3e-4)
%   'Trials'           total hop attempts                            (default 800)
%   'SigmaFnPct'       [start end] std dev for fn proposal (rel)     (default [0.10 0.01])
%   'SigmaZetaPct'     [start end] std dev for zeta proposal (rel)   (default [0.25 0.03])
%   'TemperatureScale' SA temperature scale (fraction of err)        (default 0.02)
%   'FitBaseline'      true to allow nonzero baseline (rare)         (default false)
%   'LocalRefine'      run a small LM polish each accepted step      (default true)
%
%   Peak-weighting (optional; all ignored if WeightPeaks=false):
%   'WeightPeaks'      logical, emphasize bands around initial peaks (default false)
%   'PeakBand'         [lo hi] factors around each fn0 (e.g., [0.85 1.15]) (default [0.85 1.15])
%   'PeakWeight'       weight multiplier inside those bands          (default 3.0)
%
% Output:
%   fn, zeta, A  : refined poles and residues
%   Hhat         : FRF fit at (fn,zeta,A)
%   report       : struct with rmse_before, rmse_after, accept_ratio, traj, etc.

% -------------------- parse options --------------------
p = inputParser;
p.addParameter('FRF','', @(s)ischar(s)||isstring(s));
p.addParameter('Model','auto', @(s)ischar(s)||isstring(s));
p.addParameter('DCGain',[], @(x)isnumeric(x)&&isscalar(x));
p.addParameter('Delay',0, @(x)isnumeric(x)&&isscalar(x));
p.addParameter('Areal',false, @(b)islogical(b)&&isscalar(b));
p.addParameter('MinSepPct',3e-4, @(x)isnumeric(x)&&isscalar(x)&&x>=0);
p.addParameter('Trials',800, @(x)isnumeric(x)&&isscalar(x)&&x>=1);
p.addParameter('SigmaFnPct',[0.10 0.01], @(v)isnumeric(v)&&numel(v)==2);
p.addParameter('SigmaZetaPct',[0.25 0.03], @(v)isnumeric(v)&&numel(v)==2);
p.addParameter('TemperatureScale',0.02, @(x)isnumeric(x)&&isscalar(x)&&x>0);
p.addParameter('FitBaseline',false, @(b)islogical(b)&&isscalar(b));
p.addParameter('LocalRefine',true, @(b)islogical(b)&&isscalar(b));
% Optional peak-weighting
p.addParameter('WeightPeaks',false, @(b)islogical(b)&&isscalar(b));
p.addParameter('PeakBand',[0.85 1.15], @(v)isnumeric(v)&&numel(v)==2&&all(v>0));
p.addParameter('PeakWeight',3.0, @(x)isnumeric(x)&&isscalar(x)&&x>=1);

p.parse(varargin{:});
opt = p.Results;

% -------------------- setup --------------------
f  = f(:);                 % N×1
Hc = Hc(:);                % N×1
Nm = numel(fn0);

frf    = char(opt.FRF);
model  = char(opt.Model);
tau    = opt.Delay;
dc     = opt.DCGain;
if isempty(dc)
    if strcmp(frf,''), dc = 1+0i; else, dc = 0; end
end

% prefer project helper if available
if exist('transfer_from_params','file') == 2
    xfer = @(A,fnn,zz) transfer_from_params(f, A, fnn, zz, ...
               'FRF', frf, 'Model', 'auto', 'DCGain', dc, 'Delay', tau);
else
    xfer = @(A,fnn,zz) transfer_from_params_local(f, A, fnn, zz, ...
               'FRF', frf, 'Model', 'auto', 'DCGain', dc, 'Delay', tau);
end

% initial LS residues if A0 is empty/mismatched
if isempty(A0) || numel(A0)~=Nm
    A0 = ls_residues(f, Hc, fn0, zeta0, frf, model, dc, tau, opt.Areal);
end

% --- build static weights (optional peak emphasis) ---
weights = ones(size(f));
if opt.WeightPeaks && ~isempty(fn0)
    pb = opt.PeakBand;  wpk = opt.PeakWeight;
    for j=1:numel(fn0)
        band = (f >= pb(1)*fn0(j)) & (f <= pb(2)*fn0(j));
        weights(band) = max(weights(band), wpk);
    end
end
wrmse = @(r) sqrt( mean( (weights .* abs(r)).^2 ) );

% --- helpers ---
anneal_T = @(t) opt.TemperatureScale * max(1e-6, 1 - (t-1)/max(1,opt.Trials));
lerp      = @(a,b,t) a + (b-a).*t;

% min-sep enforcement
    function [fn1, z1] = enforce_minsep(fn1, z1)
        [fn1, ord] = sort(fn1(:));
        z1 = z1(ord);
        if NumModes()>1
            minsep = max(1e-9, opt.MinSepPct) * max(1, median(fn1));
            for k=2:NumModes()
                gap = fn1(k) - fn1(k-1);
                if gap < minsep
                    mid = 0.5*(fn1(k)+fn1(k-1));
                    fn1(k-1) = mid - 0.5*minsep;
                    fn1(k)   = mid + 0.5*minsep;
                end
            end
        end
        % keep zetas in reasonable range
        z1 = max(1e-4, min(0.5, z1));
    end

    function n = NumModes(), n = numel(fn0); end

% starting state
fn   = fn0(:);
zeta = zeta0(:);
A    = A0(:);

Hhat = xfer(A, fn, zeta);
err  = wrmse(Hc - Hhat);

best.fn   = fn;  best.zeta = zeta;  best.A = A;
best.err  = err; best.Hhat = Hhat;

traj.err = zeros(opt.Trials,1);
traj.acc = false(opt.Trials,1);

% -------------------- main loop --------------------
for t = 1:opt.Trials
    % proposal scales (decay from start->end)
    sfn   = lerp(opt.SigmaFnPct(1),   opt.SigmaFnPct(2),   (t-1)/max(1,opt.Trials-1));
    szeta = lerp(opt.SigmaZetaPct(1), opt.SigmaZetaPct(2), (t-1)/max(1,opt.Trials-1));

    % propose new poles around the *current* state (SA)
    fn_p   = fn   .* (1 + sfn   * randn(Nm,1));
    zeta_p = zeta .* (1 + szeta * randn(Nm,1));
    [fn_p, zeta_p] = enforce_minsep(fn_p, zeta_p);

    % LS update residues for proposal
    A_p = ls_residues(f, Hc, fn_p, zeta_p, frf, model, dc, tau, opt.Areal);

    % model error
    H_p   = xfer(A_p, fn_p, zeta_p);
    err_p = wrmse(Hc - H_p);

    % local polish (optional)
    if opt.LocalRefine && exist('refine_modes_levenberg','file') == 2
        try
            [fn_l,z_l,A_l,H_l,rep_l] = refine_modes_levenberg( ...
                f, Hc, fn_p, zeta_p, A_p, 'FRF', frf, 'Model', model, ...
                'DCGain', dc, 'Delay', tau, 'Areal', opt.Areal, ...
                'MaxIters', 15, 'Lambda0', 1e-2, 'MinSepPct', opt.MinSepPct, ...
                'Verbose', false, 'FitBaseline', opt.FitBaseline);
            if isfield(rep_l,'rmse_after') && isfinite(rep_l.rmse_after) && rep_l.rmse_after < err_p
                fn_p=fn_l; zeta_p=z_l; A_p=A_l; H_p=H_l; err_p=rep_l.rmse_after;
            end
        catch %#ok<CTCH>
        end
    end

    % SA accept rule
    dE = err_p - err;
    T  = anneal_T(t);
    accept = (dE < 0) || (rand < exp(-dE/max(1e-12, T*max(err,1e-12))));
    if accept
        fn=fn_p; zeta=zeta_p; A=A_p; Hhat=H_p; err=err_p;
        traj.acc(t) = true;
        if err < best.err
            best.fn=fn; best.zeta=zeta; best.A=A; best.Hhat=Hhat; best.err=err;
        end
    end
    traj.err(t) = err;
end

% -------------------- results --------------------
fn   = best.fn;    zeta = best.zeta;    A = best.A;
Hhat = best.Hhat;

report = struct();
report.rmse_before = wrmse(Hc - xfer(A0(:), fn0(:), zeta0(:)));
report.rmse_after  = best.err;
report.accept_ratio = mean(traj.acc);
report.traj_err = traj.err;
report.Ntrials = opt.Trials;

end

% ======================================================================
% local helpers (used only if your project versions are absent)
function A = ls_residues(f, H, fn, zeta, frf, model, Gdc, tau, Areal)
w = 2*pi*f(:);
Nm = numel(fn);
Hcols = zeros(numel(w), Nm);
for k = 1:Nm
    Hcols(:,k) = frf_kernel_local(w, 2*pi*fn(k), zeta(k), lower(frf));
end
y = H(:);
if tau~=0, y = y .* exp(1i*w*tau); end

useMult = strcmpi(model,'multiplicative') || isempty(frf);
if useMult
    rhs = y./Gdc - 1;
else
    rhs = y - Gdc;
end

A = Hcols \ rhs;
if Areal, A = real(A); end
end

function Hhat = transfer_from_params_local(f, A, fn, zeta, varargin)
p = inputParser;
p.addParameter('DCGain',1+0i);
p.addParameter('Delay',0);
p.addParameter('FRF','');
p.addParameter('Model','auto');
p.parse(varargin{:});
G=p.Results.DCGain; tau=p.Results.Delay;
frfType = lower(char(p.Results.FRF));
modelForm = char(p.Results.Model);
if strcmpi(modelForm,'auto')
    if strcmp(frfType,''), modelForm='multiplicative'; else, modelForm='additive'; end
end

f = f(:); w = 2*pi*f;
A=A(:); fn=fn(:); zeta=zeta(:);
Nm = numel(fn);
Hcols = zeros(numel(w), Nm);
for k=1:Nm
    Hcols(:,k) = frf_kernel_local(w, 2*pi*fn(k), zeta(k), frfType);
end
if isempty(A), Hbody = zeros(size(f)); else, Hbody = Hcols*A; end
if strcmpi(modelForm,'multiplicative'), Hhat = G * (1 + Hbody);
else,                                     Hhat = G + Hbody;
end
if tau~=0, Hhat = exp(-1i*w*tau).*Hhat; end
end

function Hk = frf_kernel_local(w, wn, zeta, frfType)
den = (wn.^2 - w.^2) + 1i*(2*zeta*wn).*w;
switch frfType
    case 'disp/force', Hk = 1 ./ den;
    case 'vel/force',  Hk = (1i*w) ./ den;
    case 'acc/force',  Hk = (-w.^2) ./ den;
    otherwise,         Hk = (w.^2) ./ den;  % legacy transmissibility
end
end
