%單純形算法
%************************************************************************%
function danchunxing(func,x)
%input like this:x=[0,0;1,0;0,1],danchunxing(x);
x0=x;
u=1;r=2;e=0.5;t0=1e-10;
count=0;
xmid=zeros(1,size(x0,2));
xref=zeros(1,size(x0,2));
xdec=zeros(1,size(x0,2));
xt=zeros(1,size(x0,2));
A=zeros(1,size(x0,1));
for i=1:length(A)
    A(:,i)=func(x0(i,:));
end
x0=s_sort(func,x0,A);
xmid=fun_mid(x0);
while fun_judge(func,x0,xmid)>t0                        %判斷是否達到要求的精度
    x0=s_sort(func,x0,A);  
    xmid=fun_mid(x0);                              %求解形心
    xref=fun_ref(xmid,x0(1,:),u);                  %求解反射點
    if func(xref)<func(x0(end,:))
        xexp=fun_exp(xmid,xref,r);                 %擴張算子
        if func(xexp)<func(xref)
            x0(1,:)=xexp;
        else
            x0(1,:)=xref;
        end
    else if func(xref)<=func(x0(2,:))
            x0(1,:)=xref;
        else if func(x0(1,:))>func(xref)
                xt=xref;
            else
                xt=x0(1,:); 
            end
            xdec=fun_dec(xmid,xt,e);               %壓縮算子
            if func(xdec)<=func(xt)
                x0(1,:)=xdec;
            else
                x0=fun_x_dec(x0,x0(end,:));        %整體收縮
            end
        end
    end
    count=count+1;
end
xmin=x0(end,:),fmin=func(x0(end,:)),count

%*************************************************************************%
function x_judge=fun_judge(func,x,x_mid)                %迭代停止條件函數
x_d=0;
for i=1:size(x,1)
    x_d=x_d+(func(x(i,:))-func(x_mid))^2;
end
x_judge=sqrt(1./(size(x,2)+1)*x_d);
%*************************************************************************%
function x_s=s_sort(func,x,D)                            %數組排序函數，求解最大值
xm=zeros(1,size(x,2));                                    %最小值及次大值
for i=1:length(D)
    D(:,i)=func(x(i,:));
end
for i=1:length(D)
    for j=1:length(D)-i
        if D(j)<D(j+1)
            t=D(j);
            D(j)=D(j+1);
            D(j+1)=t;
            xm(1,:)=x(j,:);
            x(j,:)=x(j+1,:);
            x(j+1,:)=xm(1,:);
        end
    end
end
x_s=x;
%*************************************************************************%
function x_mid=fun_mid(x)                         %求解形心函數
xmt=zeros(1,size(x,2));
for i=2:size(x,1)
    xmt=xmt+x(i,:);
end
x_mid=xmt./(size(x,1)-1);
%*************************************************************************%
function x_ref=fun_ref(x_mid,x_max,u_r)           %求解反射點
x_ref=x_mid+u_r*(x_mid-x_max);
%*************************************************************************%
function x_exp=fun_exp(x_mid,x_ref,r_e)           %擴展函數
x_exp=x_mid+r_e*(x_ref-x_mid);
%*************************************************************************%
function x_dec=fun_dec(x_mid,xt,e_d)              %壓縮函數
x_dec=x_mid+e_d*(xt-x_mid);
%*************************************************************************%
function x_dec_d=fun_x_dec(x,x_min)               %收縮函數
for i=1:size(x,1)
    x(i,:)=x(i,:)+0.5*(x_min-x(i,:));
end
x_dec_d=x;
%*************************************************************************%