function [Efrac, Eabs, labels] = wavelet_component_energy_db3_parallel(X, nlevels, metric, doPlot)
% WAVELET_COMPONENT_ENERGY_DB3_PARALLEL
%   Parallel MODWT (db3) decomposition for Neuron×Time signals with per-band energy
%   and per-level histograms of the across-neuron distribution.
%
% Inputs
%   X        : (N × T) matrix, N neurons/signals by T time points
%   nlevels  : number of levels (default = 16)
%   metric   : 'energy' (default), 'l1', or 'var'  -- how to quantify "presence"
%   doPlot   : logical (default = true) to show histograms (one per level)
%
% Outputs
%   Efrac  : (N × (nlevels+1)) relative band measure per neuron
%            Columns are [D1 ... Dnlevels, A{nlevels}]. Rows ≈ sum to 1.
%   Eabs   : (N × (nlevels+1)) absolute band measure (sum of squares / L1 / var)
%   labels : 1×(nlevels+1) cell array of column labels
%
% Notes
%   - Uses MODWT + MRA (shift-invariant): modwtmra partitions the signal into
%     J detail bands + 1 approximation band at level J.
%   - Signals are mean-centered and NaNs are set to 0 prior to decomposition.
%   - Parallelization requires Parallel Computing Toolbox.
%
% Example
%   [Efrac, Eabs, labels] = wavelet_component_energy_db3_parallel(activity_Fcorr, 16);
%   % Efrac is N×17 with fractions per band; histograms are displayed automatically.

    if nargin < 2 || isempty(nlevels), nlevels = 16; end
    if nargin < 3 || isempty(metric),   metric   = 'energy'; end
    if nargin < 4 || isempty(doPlot),   doPlot   = true;     end

    validatestring(metric, {'energy','l1','var'}, mfilename, 'metric');
    assert(ismatrix(X) && isnumeric(X), 'X must be a numeric 2D matrix (Neurons × Time).');

    [N, T] = size(X);
    waveletName = 'db3';

    % Output containers
    Eabs  = zeros(N, nlevels+1);
    Efrac = zeros(N, nlevels+1);

    % Labels: D1..DJ, A{J}
    labels = [arrayfun(@(j) sprintf('D%d', j), 1:nlevels, 'UniformOutput', false), ...
              {sprintf('A%d', nlevels)}];

    % Try to start a pool (safe if already running)
    try
        if license('test','Distrib_Computing_Toolbox') %#ok<DCPT>
            p = gcp('nocreate');
            if isempty(p)
                parpool; %#ok<*NOPRT>
            end
        end
    catch ME
        warning('Parallel pool not started: %s. Falling back to serial loop.', ME.message);
    end

    % Parallel loop over neurons
    parfor i = 1:N
        x = X(i, :);

        % Mean-center; handle NaNs
        m = mean(x(~isnan(x)));
        if isnan(m), m = 0; end
        x = x - m;
        x(isnan(x)) = 0;

        % MODWT and MRA
        wt  = modwt(x, waveletName, nlevels);   % (nlevels+1) × T
        mra = modwtmra(wt, waveletName);        % (nlevels+1) × T

        % Band "presence" per level
        switch metric
            case 'energy' % sum of squares (classic energy)
                band = sum(mra.^2, 2);                  % (nlevels+1)×1
                denom = sum(x.^2);
                if denom <= eps, denom = sum(band) + eps; end  % fallback for flat signals
            case 'l1'     % total absolute amplitude
                band = sum(abs(mra), 2);
                denom = sum(band) + eps;
            case 'var'    % variance per band
                % unbiased=false (flag 0) to match energy notion of power
                band = var(mra, 0, 2) * T;              % scale like energy
                denom = sum(band) + eps;
        end

        Eabs(i, :)  = band.';
        Efrac(i, :) = (band / denom).';
    end

    % Optional plots: per-level histograms across neurons
    if doPlot
        figure('Name','Wavelet band distributions','Color','w');
        % Try to arrange in a nice grid: 17 bands -> 4×5 works (last empty), or auto
        nBands = nlevels + 1;
        nCols  = 5;
        nRows  = ceil(nBands / nCols);
        tiledlayout(nRows, nCols, 'TileSpacing','compact', 'Padding','compact');

        isFraction = strcmp(metric,'energy'); % fractions are in [0,1] only for 'energy' w/ denom=sum(x.^2)
        for k = 1:nBands
            nexttile;
            histogram(Efrac(:,k), 'NumBins', 30, 'Normalization', 'pdf');
            title(labels{k}, 'Interpreter','none');
            if isFraction, xlim([0 1]); end
            xlabel('Value'); ylabel('PDF');
            grid on;
        end
        sgtitle(sprintf('Across-neuron distribution per band (%s, db3, %d levels)', metric, nlevels));
    end
end
