function returnValue=BPSO()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% default %%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
choice=4;
kd=0.5;maxStep=2000;C1=2;C2=2;PSONUM=100;Dn=100;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% init %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%mainMatrix=[Xi,Vi,XiBest,currentFitness,LocalFitnessBest];
% Xi=[1,D],Vi=[D+1,2*D],XiBest=[2*D+1,3*D],currentFitness=3*D+1
% LobalFitnessBest=3*D+2
if(choice==1) margin=600;
elseif(choice==2) margin=30;
elseif(choice==3) margin=100;
elseif(choice==4) margin=5.12;
else margin=100;%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
end
BPSO=zeros(PSONUM,3*Dn+2);
tempMatrixCriewank=zeros(1,Dn);
tempMatrixRosenbrock=zeros(1,Dn-1);
BPSO(:,1:2*Dn)=margin*rands(PSONUM,2*Dn);
BPSO(:,Dn+1:2*Dn)=kd*BPSO(:,Dn+1:2*Dn);
BPSO(:,2*Dn+1:3*Dn)=BPSO(:,1:Dn);
if(choice==1)
    for j=1:PSONUM
        for i=1:Dn
            tempMatrixCriewank(1,i)=cos(BPSO(j,i)/sqrt(i));
        end
        BPSO(j,3*Dn+1)=sum(BPSO(j,1:Dn).^2,2)/4000-prod(tempMatrixCriewank,2)+1;
        BPSO(j,3*Dn+2)=BPSO(j,3*Dn+1);
    end
elseif(choice==2)
    BPSO(:,3*Dn+1)=exp(1)+20-20*exp(-0.2*sqrt(mean(BPSO(:,1:Dn).^2,2)))-exp(mean(cos(BPSO(:,1:Dn)*2*pi),2));
    BPSO(:,3*Dn+2)=BPSO(:,3*Dn+1);
elseif(choice==3)
    for j=1:PSONUM
        for i=1:Dn-1
            tempMatrixRosenbrock(1,i)=100*(BPSO(j,i+1)-BPSO(j,i).^2).^2+(BPSO(j,i)-1).^2;
        end
        BPSO(j,3*Dn+1)=sum(tempMatrixRosenbrock,2);
        BPSO(j,3*Dn+2)=BPSO(j,3*Dn+1);
    end
elseif(choice==4)%^^^^^^^^^^^^^
    BPSO(:,3*Dn+1)=sum(BPSO(:,1:Dn).^2+10-10*cos(BPSO(:,1:Dn)*2*pi),2);
    BPSO(:,3*Dn+2)=BPSO(:,3*Dn+1);
else%^^^^^^^^^^^^^^^^^^^^^^^^^^^
    BPSO(:,3*Dn+1)=0.5-((sin(sqrt(BPSO(:,1).^2+BPSO(:,2).^2))).^2-0.5)./(1+0.001*(BPSO(:,1).^2+BPSO(:,2).^2)).^2;
    BPSO(:,3*Dn+2)=BPSO(:,3*Dn+1);
end
[GlobalBest,GlobalIndex]=min(BPSO(:,3*Dn+2),[],1);
counter=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:maxStep
    for j=1:PSONUM
        %%%%%%%%%%%%%%%%%%%% Update Vi
        for k=1:Dn
            BPSO(j,Dn+k)=(0.9-i/maxStep*0.5)*BPSO(j,Dn+k)+rand*C1*(BPSO(j,2*Dn+k)-BPSO(j,k))+rand*C2*(BPSO(GlobalIndex,2*Dn+k)-BPSO(j,k));
            if(BPSO(j,Dn+k)>kd*margin)
                BPSO(j,Dn+k)=kd*margin;
            elseif(BPSO(j,Dn+k)<-kd*margin)
                BPSO(j,Dn+k)=-kd*margin;
            end
        end
        %%%%%%%%%%%%%%%%%%%% Update tempXi
        BPSO(j,1:Dn)=BPSO(j,1:Dn)+BPSO(j,Dn+1:2*Dn);
        for k=1:Dn
            if(BPSO(j,k)>margin)
                BPSO(j,k)=margin;
            elseif(BPSO(j,k)<-margin)
                BPSO(j,k)=-margin;
            end
        end
        %%%%%%%%%%%%%%%%%%%% Update Fitness
        if(choice==1)
            for k=1:Dn
                tempMatrixCriewank(1,k)=cos(BPSO(j,k)/sqrt(k));
            end
            BPSO(j,3*Dn+1)=sum(BPSO(j,1:Dn).^2,2)/4000-prod(tempMatrixCriewank,2)+1;
        elseif(choice==2)
            BPSO(j,3*Dn+1)=exp(1)+20-20*exp(-0.2*sqrt(mean(BPSO(j,1:Dn).^2,2)))-exp(mean(cos(BPSO(j,1:Dn)*2*pi),2));    
        elseif(choice==3)
            for k=1:Dn-1
                tempMatrixRosenbrock(1,k)=100*(BPSO(j,k+1)-BPSO(j,k).^2).^2+(BPSO(j,k)-1).^2;
            end
            BPSO(j,3*Dn+1)=sum(tempMatrixRosenbrock,2);
        elseif(choice==4)
            BPSO(j,3*Dn+1)=sum(BPSO(j,1:Dn).^2+10-10*cos(BPSO(j,1:Dn)*2*pi),2);
        else%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
            BPSO(j,3*Dn+1)=0.5-(sin(sqrt(BPSO(j,1).^2+BPSO(j,2).^2)).^2-0.5)./(1+0.001*(BPSO(j,1).^2+BPSO(j,2).^2)).^2;
        end
        %%%%%%%%%%%%%%%%%%% Update LocalBestFitness and XiBest
        if(BPSO(j,3*Dn+1)<BPSO(j,3*Dn+2))
            BPSO(j,3*Dn+2)=BPSO(j,3*Dn+1);
            BPSO(j,2*Dn+1:3*Dn)=BPSO(j,1:Dn);
        end
    end
    %%%%%%%%%%%%%%%%%%% Update GlobalBest and its Index 
    [GlobalBest,GlobalIndex]=min(BPSO(:,3*Dn+2),[],1);
    globalBestRecordMatrix(i)=GlobalBest;
end
returnValue=globalBestRecordMatrix;
% GlobalBest
% xx=1:maxStep;
% plot(xx,globalBestRecordMatrix);