function [x0_mao,gbo,gbv]=pso(x0,preparam,popparam)
%************************************
%preparam=m 預測步數
%popparam=[popsize,dimsize]
%************************************
%Initialization of PSO parameters
wmax=0.9;
wmin=0.4;
%Maximum iteration number
itmax=1000; 
c1=2;
c2=2;
m=preparam;
for iter=1:itmax
w(iter)=wmax-((wmax-wmin)/itmax)*iter;
end
%************************************
popsize=popparam(1);
dimsize=popparam(2);
[X(:,:,1),V(:,:,1),S(:,1),OUT(:,:,1)]=initpop(popsize,dimsize,x0);
% pbv=pbestvalue,pbi=pbestiter,pbo=pbestobj,gbv=gbestvalue,gbid=gbestindex,gbit=gbestiter,gbo=gbestobj
gbv=X(1,1:dimsize,1)';
gbid=1;
gbit=1;
N=1;
gbo=S(1,1);
sta_N(:,N)=[gbo,1]';
for i=1:popsize
    pbv(i,:)=X(i,1:dimsize,1);
    pbi(i)=1;
    pbo(i)=S(i,1);
    if S(i,1)<gbo
        N=N+1;
        gbo=S(i,1);
        sta_N(:,N)=[gbo,1]';
        gbv=X(i,1:dimsize,1)';
        gbid=i;
    end
end
%**************************************
%循環迭代開始，直到迭代步數結束或達到誤差閾值
for i=2:itmax
%更新位置與速度
    for j=1:popsize
        for d=1:dimsize
            V(j,d,i)=w(i)*V(j,d,i-1)+c1*rand()*(pbv(j,d)-X(j,d,i-1))+c2*rand()*(gbv(d)-X(j,d,i-1));
            X(j,d,i)=X(j,d,i-1)+V(j,d,i);
            if X(j,2,i)>1
                X(j,2,i)=1;
            else if X(j,2,i)<0
                    X(j,2,i)=0;
                end
            end
            if X(j,4,i)>1
                X(j,4,i)=1;
            else if X(j,4,i)<0
                    X(j,4,i)=0;
                end
            end
        end
    end
 %更新個體與全局最優值
    [S(:,i),OUT(:,:,i)]=objfunc(X(:,:,i),x0);
    for j=1:popsize
        if S(j,i)<pbo(j)
            pbo(j)=S(j,i);
            pbi(j)=i;
            pbv(j,:)=X(j,:,i);
        end
        if S(j,i)<gbo
            N=N+1;
            gbo=S(j,i);
            sta_N(:,N)=[gbo,i]';
            gbv=X(j,:,i)';
            gbid=j;
            gbit=i;
        end
    end
end
%*************************************
%輸出最優結果，包括預測
%Best iter=gbit;
%Best estimation error=gbo;
%Best particle=gbid;
%Best parameters producing the best prediction parameters=gbv
%Best parameters to be used predicting the expecting result=OUT(gbid,:,gbit);
%Plot gbit's error i.e. S(gbid,1:gbit) v.s. iters and
%the prediction value caculated using OUT(gbid,:,gbit) v.s. time
%caculate prediction value first,
for k=1:m
    T=OUT(gbid,:,gbit);
    %[a,u,a_xing,u_xing]=T
    a=T(1);
    u=T(2);
    a_xing=T(3);
    u_xing=T(4);
    x_0_pie(1)=x0(1)+gbv(1);
    x_1_mao(k)=(x_0_pie(1)-u/a)*exp(-a*(k-1))+u/a;
    ipsilon_0(1)=x0(1)+gbv(1)-x_1_mao(1);
    ipsilon_0_pie(1)=ipsilon_0(1)+gbv(3);
    ipsilon_0_mao(k)=(ipsilon_0_pie(1)-u_xing/a_xing)*(1-exp(a_xing))*(exp(-a_xing*(k-1)))-gbv(3);
    x_1_xing_mao(k)=x_1_mao(k)+ipsilon_0_mao(k);    
 %   x_0_pie(k)=x0(k)+gbv(1);
 %   x_1(k)=sum(x_0_pie(k));
 %   x_1_mao(k)=(x_0_pie(1)-u/a)*exp(-a*(k-1))+u/a;
  %  x_0_pie_mao(k)=(x_0_pie(1)-u/a)*(1-exp(a))*exp(-a*(k-1));
 %   x_0_mao(k)=x_0_pie_mao(k)-gbv(1);
 %   e_0(k)=x_1(k)-x_1_mao(k);
 %   e_0_pie(k)=e_0(k)+gbv(3);
  %  e_1(k)=sum(e_0_pie(k));
 %   e_1_mao(k)=(e_0_pie(1)-u_xing/a_xing)*exp(-a_xing*(k-1))+u_xing/a_xing;
 %   e_0_pie_mao(k)=(e_0_pie(1)-u_xing/a_xing)*(1-exp(a_xing))*exp(-a_xing*(k-1));
 %   e_0_mao(k)=e_0_pie_mao(k)-gbv(3);
  %  x_1_xing_mao(k)=x_1_mao(k)+e_0_mao(k);   
end
x_0_xing_mao(1)=x0(1);
for k=2:m
    x_0_xing_mao(k)=x_1_xing_mao(k)-x_1_xing_mao(k-1)-gbv(1);
end
x0_mao=x_0_xing_mao;
plot(20:N,sta_N(1,20:N),'b-')
figure(2)
plot(1:N,sta_N(2,1:N),'go');
figure(3);
plot(1:size(x0,1),x0,'b-o',1:m,x_0_xing_mao,'g-*');
%figure 4 plot several particles error behavior,chosen randomly
figure(4);
for i=1:6
    subplot(6,1,i);
    plot(1:itmax,S(i,1:itmax));
end
%output more detailed information
%'error of the last iteration:'
%gbo
%'the best results to caculate the corespondant parameters:'
%gbv(1) 
%gbv(3)
%gbv(2) 
%gbv(4)
%********************************************
%accurate test
%平均誤差
%原點誤差
%誤差和（殘差和）
%殘差平方和
%后驗差比值
%小誤差概率
n=size(x0,1);
e=abs(x0-x0_mao(1:n)')
e_re=e./x0
e_ba=sum(abs(e))/n
e_proto=abs(e(n))/x0(n)
s2=sqrt(sum(e.^2)/n)
x_ba=sum(x0)/n;
s1=sqrt(sum((x0-x_ba).^2)/n);
C=s2/s1
count=abs(e-e_ba)./s1;
k=0;
for i=1:n
    if count(i)<0.6745
        k=k+1;
    end
end
p=k/n