function [w,q,best] = DEGoStrain(Parameter,net,train_data,train_target)
%% parameter input
F = Parameter.F;
CR = Parameter.CR;
popsize = Parameter.popsize;
D = Parameter.D;
iter = 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
% fitness = inf(popsize,1);
convergence = [];
%%
for i = 1:iter
%     if i>6 && std(convergence((i-6):(i-1)))==0
%         re_de = 0;
%     else
%         re_de = 1;
%     end
    % evalution function
    U = pop;
    fit_U = fitness;
    [~,k] = min(fit_U);
    [r1, r2, r3] = getindex(popsize);
    % Implement DE/rand/1 mutation
    V = pop(r1, :) + F * (pop(r2, :) - pop(r3, :));
    
    % Check whether the mutant vector violates the boundaries or not
    V = repair(V, lu(1,1),lu(2,1));
    for j = 1:popsize
%         if  re_de == 0 && j~=k
%             radius = V(j,:) - pop(k,:);%run away from the population best
%             U(j,:) = V(j,:) + (rand(1,1)-0.5) * radius;    %Chaos_p(c_j,c_d)
%         else
            j_rand = floor(rand * D) + 1;
            t = rand(1, D) < CR;
            t(1, j_rand) = 1;
            t_ = 1 - t;
            U(j, :) = t .* V(j, :) + t_ .* pop(j, :);
%         end
    end
    
    %     fitness = evolution(Parameter,U,train_data,train_target,net);
    
    % 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);
    end
    
    % Select the better one between the target vector and its trial vector
    for j = 1:popsize
        if fit_U(j, :) <= fitness(j, :)
            pop(j, :) = U(j, :);
            fitness(j, :) = fit_U(j, :);
        end
    end
    best = min(fitness);
    convergence= [convergence; best];    
    %     disp(['The iteration ' num2str(i) ' error:' num2str(best)]);
end
%%
[~,index] = min(fitness);
best_population = pop(index,:);
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;
end