function [w,q,best,out] = bbotrain(Parameter,net,train_data,train_target)
%% parameter input
F = Parameter.F;
CR = Parameter.CR;
popsize = Parameter.popsize;
D = Parameter.D;
MaxIt = Parameter.iter;

M = net.M;
fitness = zeros(popsize,1);
[~,J] =size(train_data); % J: the dimension of samples

%% initial population
lu = [-1 * ones(1, D); 1 * ones(1, D)];
pop = repmat(lu(1, :), popsize, 1) + rand(popsize, D) .* (repmat(lu(2, :) - lu(1, :), popsize, 1));

% objective function

% fitness = evolution(Parameter,pop,train_data,train_target,net);
for popindex = 1:popsize
    w = pop(popindex,1:J*M);
    q = pop(popindex,J*M+1:D);
    w=reshape(w,J,M);
    q=reshape(q,J,M);
    net.w = w;
    net.q = q;
    train_fit = my_DNM(train_data,net);
    cost = (train_fit - train_target).^2;
    fitness(popindex) = mean(cost);
end

%%
%bbo train
KeepRate=0.2;     % Keep Rate
VarSize=[1 D];
VarMin = lu(1,1);
VarMax = lu(2,1);
nKeep=round(KeepRate*popsize);     % Number of Kept Habitats

nNew=popsize-nKeep;                % Number of New Habitats

% Migration Rates
mu=linspace(1,0,popsize);          % Emmigration Rates
lambda=1-mu;                    % Immigration Rates

alpha=0.9;

pMutation=0.1;

sigma=0.02*(lu(2,1)-lu(1,1));

habitat.Position=[];
habitat.Cost=[];

% Create Habitats Array
pop=repmat(habitat,popsize,1);

% Initialize Habitats
for ii=1:popsize
    pop(ii).Position=unifrnd(VarMin,VarMax,VarSize);
    U(ii,:) = pop(ii).Position;
end

% objective function
for popindex = 1:popsize
    w = U(popindex,1:J*M);
    q = U(popindex,J*M+1:D);
    w=reshape(w,J,M);
    q=reshape(q,J,M);
    net.w = w;
    net.q = q;
    train_fit = my_DNM(train_data,net);
    cost = (train_fit - train_target).^2;
    fit_U(popindex,:) = mean(cost);
    pop(popindex).Cost = fit_U(popindex,:);
end

% Sort New Population
[~, SortOrder]=sort([pop.Cost]);
pop=pop(SortOrder);

for it = 1:MaxIt
    newpop=pop;
    for i=1:popsize
        for k=1:D
            % Migration
            if rand<=lambda(i)
                % Emmigration Probabilities
                EP=mu;
                EP(i)=0;
                EP=EP/sum(EP);
                
                % Select Source Habitat
                j=RouletteWheelSelection(EP);
                
                % Migration
                newpop(i).Position(k)=pop(i).Position(k) ...
                    +alpha*(pop(j).Position(k)-pop(i).Position(k));
            end
            % Mutation
            if rand<=pMutation
                newpop(i).Position(k)=newpop(i).Position(k)+sigma*randn;
            end
        end
        % Apply Lower and Upper Bound Limits
        newpop(i).Position = max(newpop(i).Position,VarMin);
        newpop(i).Position = min(newpop(i).Position,VarMax);
        U(i,:) = newpop(i).Position;
    end
    % objective function
    
    for popindex = 1:popsize
        w = U(popindex,1:J*M);
        q = U(popindex,J*M+1:D);
        w=reshape(w,J,M);
        q=reshape(q,J,M);
        net.w = w;
        net.q = q;
        train_fit = my_DNM(train_data,net);
        cost = (train_fit - train_target).^2;
        fit_U(popindex) = mean(cost);
        newpop(popindex).Cost = fit_U(popindex);
    end
    
    % Sort New Population
    [~, SortOrder]=sort([newpop.Cost]);
    newpop=newpop(SortOrder);
    % Select Next Iteration Population
    pop=[pop(1:nKeep)
        newpop(1:nNew)];
    % Sort Population
    [~, SortOrder]=sort([pop.Cost]);
    pop=pop(SortOrder);
    BestCost(it)=pop(1).Cost; 
     disp(['it=',num2str(it),'|MSE=',num2str(BestCost(it))]);
end
%%

[best,index] = min(pop(1,:).Cost);
best_population = newpop(index,:).Position;
w = best_population(1:J*M);
q = best_population(J*M+1:D);
w=reshape(w,J,M);
q=reshape(q,J,M);
net.w = w;
net.q = q;
out.convergence= BestCost;
out.w=w;
out.q=q;
best = BestCost(MaxIt);
end
