function [Fcorr, alpha_vec, corr_vec, energy_preserved, info] = ...
    correct_neuropil_dc_safe(F, Fneu, baseline_ends, varargin)
% CORRECT_NEUROPIL_DC_PAR  Neuropil subtraction that preserves DC (no smoothing).
%
%   [Fcorr, alpha_vec, corr_vec, energy_preserved, info] = ...
%       correct_neuropil_dc_par(F, Fneu, baseline_ends, 'Name',Value,...)
%
% Inputs
%   F            : N×T fluorescence (neurons × time)
%   Fneu         : N×T neuropil (same size as F)
%   baseline_ends: if scalar B -> baseline = 1:B; else vector/logical indices of baseline frames
%
% Name-Value (optional)
%   'amax'        : cap on alpha (default 0.7)
%   'UseParallel' : true/false to prefer parfor (default true)
%
% Outputs
%   Fcorr            : N×T, corrected signal (DC preserved)
%   alpha_vec        : N×1, per-neuron alpha (slope)
%   corr_vec         : N×1, corr(Fneu, Fcorr) on baseline (diagnostic)
%   energy_preserved : N×1, Ecorr/Eorig using demeaned energies
%   info             : struct with base_idx, amax, usedParallel

% ---------- parse ----------
p = inputParser;
addParameter(p,'amax',0.7);
addParameter(p,'UseParallel',true);
parse(p,varargin{:});
amax   = p.Results.amax;
wantPar= p.Results.UseParallel;

% ---------- prep ----------
F    = double(F);
Fneu = double(Fneu);
[N,T]= size(F);
assert(isequal(size(Fneu),[N,T]), 'F and Fneu must have same size.');

% baseline indices
if isscalar(baseline_ends)
    base_idx = 1:max(1, round(baseline_ends));
else
    if islogical(baseline_ends), base_idx = find(baseline_ends);
    else, base_idx = baseline_ends(:)'; end
end
base_idx = base_idx(base_idx>=1 & base_idx<=T);
assert(~isempty(base_idx), 'baseline_ends defines no valid frames.');

% ---------- outputs ----------
Fcorr            = zeros(N,T);
alpha_vec        = zeros(N,1);
corr_vec         = zeros(N,1);
energy_preserved = zeros(N,1);

% ---------- loop (parfor if available) ----------
usePar = wantPar && hasParallelToolbox();
if usePar
    parfor i = 1:N
        [Fcorr(i,:), alpha_vec(i), corr_vec(i), energy_preserved(i)] = ...
            subtract_one(F(i,:), Fneu(i,:), base_idx, amax);
    end
else
    for i = 1:N
        [Fcorr(i,:), alpha_vec(i), corr_vec(i), energy_preserved(i)] = ...
            subtract_one(F(i,:), Fneu(i,:), base_idx, amax);
    end
end

info = struct('base_idx',base_idx,'amax',amax,'usedParallel',usePar);
end

% ================= Worker (per neuron) =================
function [Fc, a, rcorr, ekeep] = subtract_one(Fi, Fneui, base_idx, amax)
% --- Center neuropil on BASELINE (fit domain = raw; no smoothing)
nb_fit     = Fneui(base_idx);
mu_neu_fit = median(nb_fit,'omitnan');      % baseline center
x          = nb_fit - mu_neu_fit;          % zero-mean regressor on baseline
y          = Fi(base_idx);

% --- OLS with intercept on baseline, then cap
x = x(:); y = y(:);
m = isfinite(x) & isfinite(y);
if nnz(m) >= 3 && (max(x(m)) - min(x(m))) > 0
    X = [ones(nnz(m),1), x(m)];            % [intercept, slope]
    b = X \ y(m);
    a = b(2);                               % slope = alpha
else
    a = 0;
end
if ~isfinite(a), a = 0; end
a = max(0, min(a, amax));                   % clamp 0 ≤ α ≤ amax

% --- Build centered neuropil on RAW domain (preserve DC)
mu_neu_raw = mu_neu_fit;                    % same as fit center (no smoothing case)
nreg       = Fneui - mu_neu_raw;            % zero-mean on baseline

% --- Subtract ONLY α·(centered neuropil)  → DC preserved
Fc = Fi - a * nreg;

% --- Diagnostics
rcorr = corr(Fneui(base_idx)', Fc(base_idx)', 'rows','pairwise');

% Demeaned energies (event power) before/after
y0    = Fi - mean(Fi,'omitnan');
yc    = Fc - mean(Fc,'omitnan');
Eorig = nansum(y0.^2);
Ecorr = nansum(yc.^2);
ekeep = (Eorig>0) * (Ecorr / max(Eorig, eps));
end

% ===== helper =====
function tf = hasParallelToolbox()
    tf = license('test','Distrib_Computing_Toolbox');
end
