function [sFeat, Sf, Nf, curve] = jWOA(feat, label, N, max_Iter, HO)

% Whale Optimization Algorithm (WOA) for Feature Selection
% Inputs:
%   - feat: Feature matrix
%   - label: Class labels
%   - N: Population size
%   - max_Iter: Maximum iterations
%   - HO: Hold-out validation partition
% Outputs:
%   - sFeat: Selected features
%   - Sf: Indices of selected features
%   - Nf: Number of selected features
%   - curve: Convergence curve of best fitness

% Parameter Initialization
lb = 0; % Lower bound
ub = 1; % Upper bound
b = 1;  % Spiral movement parameter

dim = size(feat, 2); % Number of features
fun = @jFitnessFunction; % Fitness function

% Initialize search agents
X = lb + (ub - lb) * rand(N, dim);

% Initialize fitness values
fit = zeros(1, N);
fitG = inf; % Global best fitness

% Evaluate initial solutions
for i = 1:N
    fit(i) = fun(feat, label, (X(i, :) > 0.5), HO);
    if fit(i) < fitG
        fitG = fit(i); 
        Xgb = X(i, :); % Global best position
    end
end

curve = inf(1, max_Iter); % Initialize convergence curve

t = 1; % Iteration counter

% Main WOA optimization loop
while t <= max_Iter
    a = 2 - t * (2 / max_Iter); % Linearly decreasing parameter
    
    for i = 1:N
        A = 2 * a * rand() - a; % Exploration coefficient
        C = 2 * rand(); % Random coefficient
        p = rand(); % Probability factor
        l = -1 + 2 * rand(); % Random number in range [-1,1]

        if p < 0.5 % Exploration and exploitation phase
            if abs(A) < 1
                % Encircling prey (exploitation)
                for d = 1:dim
                    Dx = abs(C * Xgb(d) - X(i, d));
                    X(i, d) = Xgb(d) - A * Dx;
                end
            else
                % Searching for prey (exploration)
                k = randi([1, N]); % Select a random search agent
                for d = 1:dim
                    Dx = abs(C * X(k, d) - X(i, d));
                    X(i, d) = X(k, d) - A * Dx;
                end
            end
        else
            % Bubble-net attacking (spiral movement)
            for d = 1:dim
                dist = abs(Xgb(d) - X(i, d));
                X(i, d) = dist * exp(b * l) * cos(2 * pi * l) + Xgb(d);
            end
        end
        
        % Ensure boundaries are respected
        X(i, X(i, :) > ub) = ub;
        X(i, X(i, :) < lb) = lb;
    end
    
    % Evaluate new solutions
    for i = 1:N
        fit(i) = fun(feat, label, (X(i, :) > 0.5), HO);
        if fit(i) < fitG
            fitG = fit(i);
            Xgb = X(i, :);
        end
    end
    
    % Store best fitness in convergence curve
    curve(t) = fitG;
    t = t + 1;
end

% Extract selected features
Pos = 1:dim;
Sf = Pos((Xgb > 0.5) == 1);
Nf = length(Sf);
sFeat = feat(:, Sf); % Final selected feature subset

end