function [RBF_surrogate,app_fobj,Trainset_size] = globsurrogate_oper_vf7(mode,iter,currFEs,FEmax,old_surrogate,glob_Pos,glob_Fit,swarmPos_curr)
%% globsurrogate_oper_vf6 uses a different objective function approximation using Radial Base Functions (RBF)
% alternative to newrb by using: rbfcreate and rbfinterp functions
minPos_surr=30; % minimum allowed number of positions required to assemble the surrogate
maxPos_surr=150; % maximum allowed number of positions required to assemble the surrogate
switch mode
    %________________________________________________________________________________%
    case 0 %Build Surrogate Model approximation
        %________________________________________________________________________________%
        if (isempty(swarmPos_curr)) %situation that occurs on the 1st iteratiomn/1st generation swarm (recently generated positions/search agents)
            x_conf=glob_Pos;
            fitness_conf=glob_Fit;
        else %from the second generation onward it is likely that onlyt part of the data in the global database will be employed
            %->The global surrogate model is bulit using the data samples in the whole decision space in the beggining and then the
            %search proceeds only part of the samples in the global database, that, are near the location of the current swarm, to train the model
            %-------------------------------------------------------------------------%
            alpha_coef=0.01+0.305*exp(-0.0015*iter); %alpha variable that determines how far beyond the current swarm the subspace finds positions in the DB (ensures that even if many iterations are made, the spread never goes below 0.05)
            %-------------------------------------------------------------------------%
            % find the current lower and upper bounds of the current swarm and the global database (all evaluated positions with the real fobj
            xbounds_d_curr=[min(swarmPos_curr,[],1); max(swarmPos_curr,[],1)]; % minimum and maximum values of the current swarm on the d.th dimension: min_d_curr (1st line) and max_d_curr (2nd line)
            xbounds_d_all=[min(glob_Pos,[],1); max(glob_Pos,[],1)];  % minimum and maximum values of the whole search space (all positions on the global data) on the d.th dimension: min_d_all (1st line) and max_d_all (2nd line)
            %-------------------------------------------------------------------------%
            % define the searchspace that builds the global surrogate outside the space defined by the current swarm, and the spread coef defines how far it searches
            aux_spmin=xbounds_d_curr(1,:)-alpha_coef*(xbounds_d_curr(2,:)-xbounds_d_curr(1,:));
            aux_spmax=xbounds_d_curr(2,:)+alpha_coef*(xbounds_d_curr(2,:)-xbounds_d_curr(1,:));
            sp_xbounds_d=[max(aux_spmin,xbounds_d_all(1,:)); min(aux_spmax,xbounds_d_all(2,:))]; % subspace (bounds) of datasamples to build the global surrogate
            %-------------------------------------------------------------------------%
            % with the subspace defined, assemble (gather the positions inside this subspace) the training pair of data samples
            aux=1;
            x=[];
            fitness=[];
            sel_idxs=[];
            for i=1:size(glob_Pos,1) %rotate all the positions/search agents stored in the DB
                if (~(sum(glob_Pos(i,:)<sp_xbounds_d(1,:)) || sum(glob_Pos(i,:)>sp_xbounds_d(2,:)))) %i.e. if position inside the choosen subspace: does not violate the subspace bounds -> ~ ( (sum(x<lb) || sum(x>ub) )
                    sel_idxs=[sel_idxs, i];
                    x(aux,:)=glob_Pos(i,:);
                    fitness(aux,1)=glob_Fit(i);
                    aux=aux+1;
                end
            end
            %-------------------------------------------------------------------------%
            % check the integrity of the assembled training data
            if (isempty(x))  % if the x and consequently fitness matrix (supposedly contains the  selected positions inside the subspace) are empty
                % than select the closest (to the current swarm) "minPos_surr" positions  
                % i.e., the minimum allowable number of positions from the DB, to assemble the DB
%                 % Alternative 1
%                 mean_swarmPos_curr=mean(swarmPos_curr); %mean current swarm position
%                 dist_swarmPos_curr2DB=zeros(size(glob_Pos,1),1);
%                 for j=1:size(glob_Pos,1) %calculate the distance of each stored DB position to the mean current swarm position
%                     dist_swarmPos_curr2DB(j)=norm(glob_Pos(j,:)-mean_swarmPos_curr);
%                 end
%                 %$
                % Alternative 2
                mean_swarmPos_curr=repmat(mean(swarmPos_curr,1),size(glob_Pos,1),1); %mean current swarm position (repeated alongside the dimensions of the DB
                dist_swarmPos_curr2DB=vecnorm((glob_Pos-mean_swarmPos_curr)'); % If input is a matrix, then vecnorm returns the norm of each column, hence the '
                %$
                [~,sortdist_idxs]=sort(dist_swarmPos_curr2DB,'ascend');
                x_conf=glob_Pos(sortdist_idxs(1:minPos_surr),:);
                fitness_conf=glob_Fit(sortdist_idxs(1:minPos_surr));
            elseif (size(x,1)<minPos_surr) %or if it does not contain enough (the minimum number of positions to assemble the surrogate)
                misslength=minPos_surr-size(x,1); %number of missing positions
                x_conf=[x; ones(misslength,size(x,2))]; %add the missing ones to reach the minimum allowable number of positions (without repetition)
                fitness_conf=[fitness; ones(misslength,1)];
%                 % Alternative 1
%                 mean_swarmPos_curr=mean(swarmPos_curr); %mean current swarm position
%                 dist_swarmPos_curr2DB=zeros(size(glob_Pos,1),1);
%                 for j=1:size(glob_Pos,1) %calculate the distance of each stored DB position to the mean current swarm position
%                     dist_swarmPos_curr2DB(j)=norm(glob_Pos(j,:)-mean_swarmPos_curr);
%                 end
%                 %$
                % Alternative 2
                mean_swarmPos_curr=repmat(mean(swarmPos_curr,1),size(glob_Pos,1),1); %mean current swarm position (repeated alongside the dimensions of the DB
                dist_swarmPos_curr2DB=vecnorm((glob_Pos-mean_swarmPos_curr),2,2); % vecnorm returns the norm type 2 of each row (2)
                %$
                [~,sortdist_idxs]=sort(dist_swarmPos_curr2DB,'ascend');
                for j=1:length(sortdist_idxs) %check the agents stored in the DB and find non repeating with the 
                    % minimum distance to the mean curr swarm position ones until enoughh are gathered
                    if (~sum(sortdist_idxs(j)==sel_idxs) && (misslength>0)) %if the index sortdist_idxs(j) position was not yet included 
                        % in the selected training set (x)  
                        x_conf(minPos_surr-misslength+1,:)=glob_Pos(sortdist_idxs(j),:);
                        fitness_conf(minPos_surr-misslength+1)=glob_Fit(sortdist_idxs(j),:);
                        misslength=misslength-1;
                    elseif (~misslength) %equal to (misslength==0)
                        break; %the minimum number of (all different) positions in the train set was reached 
                    end
                end
            elseif (size(x,1)>maxPos_surr) % if the training set is too large will it will compromise the time to train, 
                % so select the top best maximum allowed number of positions
                %                     size(x,1)
                %                     idx=randperm(popsize*5);
                [~,idx]=sort(fitness,'ascend');
                x_conf=x(idx(1:maxPos_surr),:); % Alternative select the top fittest positions
                fitness_conf=fitness(idx(1:maxPos_surr));
            else % in the case that everything is ok!
                x_conf=x;
                fitness_conf=fitness;
            end
        end
        %-------------------------------------------------------------------------%
        %         length(x_conf)
        %         figure(3),plot(x_conf(:,1),x_conf(:,2),'*')
        %         pause(0.5)
        %         clf
        %________________________________________________________________________________%
        % Interpolate selected subspace data with a RBFNN to construct the global surrogate
        x_train=x_conf'; %formato cada sample/position (all dimensions) is represented by a column
        Trainset_size=length(x_train);
        y_train=fitness_conf'; %correspodente em termos de funcao objetivo
        % Radial basis networks parameters
        goal_max=0.1; %MSE goal
        goal_min=0.01; %MSE goal
        goal=goal_max-((goal_max-goal_min)/FEmax)*currFEs;
        %         spread=1; %Spread of radial basis functions. default->1
        adaptativeSPREAD=min(max((max(x_conf,[],1)-min(x_conf,[],1)),0),[],2);
        if (adaptativeSPREAD==0)
            adaptativeSPREAD=0.05;
            disp('check what happened')
        end
        spread=adaptativeSPREAD;
        %RBF interpolation: give exact value of the function at the nodes. 
        % If input data is noisy or exact values at the nodes are not desired smoothed approximation can be used instead of interpolation.
%         RBF_surrogate=rbfcreate(x_train,y_train,'RBFFunction', 'cubic','RBFConstant',spread,'RBFSmooth',goal);
        [~,RBF_surrogate]=evalc('rbfcreate(x_train,y_train,''RBFFunction'',''cubic'',''RBFConstant'',spread,''RBFSmooth'',goal)');
        % y_train: dimensions->[dim,numtrainsamples]; x_train: dimensions->[1,numtrainsamples]
        % RBF_surrogate=coeffs: Interpolation coefficients used in function rbfinterp 
        app_fobj=[];
        %________________________________________________________________________________%
    case 1 % Test/Use Approximated Surrogate Model to obtain an approximation for the objective value
        app_fobj=rbfinterp(swarmPos_curr',old_surrogate);
        RBF_surrogate=[];
        Trainset_size=[];
        %________________________________________________________________________________%
end
end
%____________________________________________________________________________________________________________________________________________________________________________________________________________________________________________%
% % Test accuracy with a new random generated particle
% newPosition=Lb+(Ub-Lb).*rand(1,dim);
% y_est=sim(net,newPosition') %estimativa usando a interpolacao com a rede
% y_real=fobj(newPosition)
% error=y_real-y_est
%____________________________________________________________%
% by Pedro Bento November 2023
% When using please cite and refer this code and its author!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%