%% ========================================================================
%  BSO OPTIMIZER FOR 23 BENCHMARK FUNCTIONS (F1–F23)
%
%  Developed by:
%   - CEATS, Ho Chi Minh City Open University, Vietnam
%   - https://ceats.ou.edu.vn/us/codes.html 
%    Link to download article: https://doi.org/10.1016/j.asoc.2026.114975 
%  Purpose:
%   - Solve a selected benchmark function (F1..F23) using a B-spline based
%     optimizer with exploration/exploitation strategies controlled by xi.
%
%  Required functions (must be in the same folder / MATLAB path):
%   - Get_Functions_details(fun)           -> [lb, ub, dim, fobj]
%   - DersOneBasisFun(p, U, i, u, n)       -> B-spline basis + derivatives
%   - boundary_condition_BSO(x, Lb, Ub, dim) -> clamp/repair to bounds
%
%  Output:
%   - G_best.pos  : best solution vector
%   - G_best.cost : best fitness value
%   - best_GridPoints_cost(k) : convergence history
%% ========================================================================

clc; clear; close all; format short e

%% -------------------- 1) SELECT BENCHMARK FUNCTION ---------------------
fun = 'F18';                                  % Benchmark name (F1–F23)
disp(['Function ', fun]);

% Get lower/upper bounds, dimension, and function handle
[lb, ub, dim, fobj] = Get_Functions_details(fun);

nvar = dim;                                   % Problem dimension
Lb   = lb .* ones(1, nvar);                   % Lower bound vector
Ub   = ub .* ones(1, nvar);                   % Upper bound vector

NP    = 30;                                   % Population size
cycle = 1000;                                 % Max iterations

%% -------------------- 2) INITIALIZE POPULATION --------------------------
% Each particle/agent has: position, velocity, learning rate (lr), cost
for q = 1:NP
    % NOTE: keeps your original initialization style
    GridPoints(q).pos  = Lb + rand(1).*(Ub - Lb);
    GridPoints(q).vel  = zeros(1, nvar);
    GridPoints(q).lr   = ones(1, nvar).*0.9;
    GridPoints(q).cost = fobj(GridPoints(q).pos);

    current_cost(q,:) = GridPoints(q).cost;   % store current fitness
    current_pos(q,:)  = GridPoints(q).pos;    % store current position
end

% Copy population for self-improvement steps
B_spline = GridPoints;

% Identify global best and worst solutions
[~, index_min] = min([GridPoints.cost]);
[~, index_max] = max([GridPoints.cost]);
G_best  = GridPoints(index_min);
G_worst = GridPoints(index_max);

% Sort (optional, used for selection steps)
[value_sort, index_sort] = sort([GridPoints.cost]); %#ok<ASGLU>

%% -------------------- 3) PRE-ALLOCATE MEMORY ----------------------------
pre_pos  = zeros(NP, nvar);                   % previous position (per agent)
pre_cost = zeros(NP, 1);                      % previous fitness (per agent)

%% -------------------- 4) DEFINE KNOT VECTOR (QUADRATIC BASIS) -----------
% Your original knot vector (degree computed later from NCtrlPts)
KntVect = [0 0 0 1 2 3 4 4 4]./4;
KntVect = KntVect ./ max(KntVect);

%% ============================ MAIN LOOP =================================
for k = 1:cycle

    % -------------------- Loop over agents ------------------------------
    for j = 1:NP

        xi = rand;                            % random switch for strategies

        % Store previous (for learning-rate adaptation)
        pre_cost(j,:) = current_cost(j,:);
        pre_pos(j,:)  = current_pos(j,:); 

        % Recompute sort (kept as in your code)
        [value_sort, index_sort] = sort([GridPoints.cost]); 

        % Random indices used to pick candidate control points
        random_indices = randperm(NP);

        % -------------------- 5) BUILD CONTROL POINT SET P_xi ------------
        % P_xi contains 6 control points: worst, some random-worst,
        % some random-best, and best 
        P_xi(1,:) = G_worst.pos;
        P_xi(2,:) = GridPoints(random_indices(end-1)).pos;
        P_xi(3,:) = GridPoints(random_indices(end-2)).pos;
        P_xi(4,:) = GridPoints(random_indices(2)).pos;
        P_xi(5,:) = GridPoints(random_indices(1)).pos;
        P_xi(6,:) = G_best.pos;

        NCtrlPts = size(P_xi, 1);

        % Degree inferred from knot vector and number of control points
        p_order = numel(KntVect) - NCtrlPts - 1;

        % -------------------- 6) DERIVATIVE AT "BEST" PARAMETER ----------
        % Compute first derivative of basis at xi=1 (as you did for "best")
        N0best = zeros(2, NCtrlPts);
        for i = 1:NCtrlPts
            N0best(:, i) = DersOneBasisFun(p_order, KntVect, i, 1, 1);
        end

        Der_B_spline_pos_best = N0best(2,:) * P_xi;          % derivative
        Der_B_spline_pos_best = boundary_condition_BSO( ...
            Der_B_spline_pos_best, Lb, Ub, dim);

        % -------------------- 7) BASIS + DERIVATIVE AT RANDOM xi ---------
        N0 = zeros(2, NCtrlPts);
        for i = 1:NCtrlPts
            N0(:, i) = DersOneBasisFun(p_order, KntVect, i, xi, 1);
        end

        % Current B-spline point and its derivative
        B_spline(j).pos        = N0(1,:) * P_xi;
        Der_B_spline(j).pos    = N0(2,:) * P_xi;
        Der_B_spline(j).pos    = boundary_condition_BSO( ...
            Der_B_spline(j).pos, Lb, Ub, dim);

     
        momentum = 0.1;

        % -------------------- 8) STRATEGY SWITCHING BY xi ----------------
        % (A) Exploitation mode 1
        if xi >= 0.75
            B_spline(j).pos = G_best.pos + (-1 + 2*rand(1,nvar)) .* ...
                ( G_best.pos - 2.*rand .* B_spline(j).lr .* Der_B_spline_pos_best );
        end

        % (B) Exploitation mode 2
        if (0.5 <= xi) && (xi < 0.75)
            B_spline(j).pos = B_spline(j).pos - (-1 + 2*rand(1)) .* ...
                B_spline(j).pos .* B_spline(j).lr .* Der_B_spline(j).pos;

            B_spline(j).pos = G_best.pos + (-1 + 2*rand(1,nvar)) .* ...
                ( G_best.pos - B_spline(j).pos );
        end

        % (C) Exploration mode 1 (velocity update + perturbation)
        if (0.25 <= xi) && (xi < 0.5)
            B_spline(j).vel = momentum .* B_spline(j).vel + (-1 + 2*rand) .* ...
                B_spline(j).lr .* Der_B_spline(j).pos;

            B_spline(j).pos = B_spline(j).pos + B_spline(j).vel;

            B_spline(j).pos = 2.*rand .* B_spline(j).pos + (-1 + 2*rand(1,nvar)) .* ...
                ( G_best.pos - B_spline(j).pos );
        end

        % (D) Exploration mode 2 (randomized mixing with bounds/best/worst)
        if xi < 0.25
            if rand > 0.5
                B_spline(j).pos = 2.*rand .* B_spline(j).pos ...
                    + 2*rand(1,nvar) .* (B_spline(j).pos - 2*rand .* G_worst.pos) ...
                    + 2*rand(1,nvar) .* (G_best.pos - 2*rand(1).*B_spline(j).pos);
            else
                B_spline(j).pos = 2.*rand .* Lb ...
                    + (B_spline(j).pos - 2.*rand(1).*Lb) ...
                    + (Ub - 2.*rand(1).*B_spline(j).pos);
            end
        end

        % -------------------- 9) BOUNDARY HANDLING + FITNESS --------------
        B_spline(j).pos  = boundary_condition_BSO(B_spline(j).pos, Lb, Ub, dim);
        B_spline(j).cost = fobj(B_spline(j).pos);

        % -------------------- 10) LEARNING-RATE ADAPTATION ----------------
        % If new trial is worse than previous cost -> decrease lr, else increase
        if B_spline(j).cost > pre_cost(j,:)
            B_spline(j).lr = B_spline(j).lr * 0.9;
        else
            B_spline(j).lr = B_spline(j).lr * 1.01;
        end

        % -------------------- 11) GREEDY REPLACEMENT ----------------------
        if B_spline(j).cost < GridPoints(j).cost
            GridPoints(j) = B_spline(j);
        end

        % Update global best / worst
        if GridPoints(j).cost < G_best.cost
            G_best = GridPoints(j);
        end
        if GridPoints(j).cost > G_worst.cost
            G_worst = GridPoints(j);
        end

        % Store for next iteration
        current_cost(j,:) = GridPoints(j).cost;
        current_pos(j,:)  = GridPoints(j).pos;

    end % end loop j (agents)

    % Save convergence history
    best_GridPoints_cost(k,:) = G_best.cost;
    best_GridPoints_pos(k,:)  = G_best.pos;

end % end loop k (iterations)

%% -------------------- 12) DISPLAY & PLOT RESULTS ------------------------
disp(['Best solution: ', num2str(G_best.pos)]);
disp(['Best fitness:  ', num2str(G_best.cost)]);

figure;
plot(1:cycle, best_GridPoints_cost, '-r', 'LineWidth', 2);
xlabel('Iteration',    'FontName', 'Times', 'FontSize', 16);
ylabel('Best fitness', 'FontName', 'Times', 'FontSize', 16);
set(gca, 'FontName', 'Times', 'FontSize', 16);
legend('BSO-optimizer', 'FontName', 'Times', 'FontSize', 16);