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

% RNFS: A robust nature-inspired feature selection for remote sensing image classification.
% 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
MinT = inf; % Initial minimum fitness value

% Fox jump parameters
c1 = 0.18; % Range [0, 0.18]
c2 = 0.82; % Range [0.19, 1]

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

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

% Evaluate initial solutions
fit = zeros(1, N);
fitG = inf;
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

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

% Main RNFS optimization loop
while iter <= max_Iter
    a = 2 * (1 - (iter / max_Iter)); % Adaptive parameter
    
    for i = 1:N
        r = rand;
        p = rand();
        p1 = rand();

        if r >= 0.5
            if p < 0.4 % Helix move
                b = 1;
                lw = -1 + 2 * rand();
                for d = 1:dim
                    dist = abs(Xgb(d) - X(i, d));
                    X(i, d) = dist * exp(b * lw) * cos(2 * pi * lw) + Xgb(d);
                end
            else % Shifts in two directions
                Time = rand(1, dim);
                sps = Xgb ./ Time;
                Distance_S_Travel = sps .* Time;
                Distance_Fox_Rat = 0.50 * Distance_S_Travel;
                tt = sum(Time) / dim;
                t = tt / 2;
                Jump = 0.5 * 9.81 * t^2;
                
                if p1 >= 0.18
                    X(i, :) = Distance_Fox_Rat .* Jump * c1;
                else
                    X(i, :) = Distance_Fox_Rat .* Jump * c2;
                end
            end
        else % Exploration phase
            X(i, :) = Xgb + randn(1, dim) .* (MinT * a);
        end

        % Boundary handling
        X(i, X(i, :) > ub) = ub;
        X(i, X(i, :) < lb) = lb;
    end
    
    % Evaluate 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
    
    curve(iter) = fitG;
    iter = iter + 1;
end

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

end