function [bestx,bestfx] = sceua(x0,bl,bu,maxn,kstop,pcento,peps,ngs,iseed,iniflg);
% Definition:
% xo ：待优化参数的初始值
% b1 : 待优化参数的下限
% bu ： 待优化参数的上限
% maxn ：进化过程中函数最大调用次数
% kstop ：集中之前最大演化数
% pcento ：允许在kstop循环前收敛的百分比变化
% iseed ：生成随机数种子
% iniflg ：初始参数矩阵的标志
% ngs ：复合型个数
% npg ：复合型顶点数目
% nps ：简单型
% nspl ：每个复合型进化的迭代数
% mings ：进化过程复合型的最小数目

global BESTX BESTF ICALL PX PF
nopt=length(x0); %点的维度
npg=2*nopt+1;%复合型顶点数目
nps=nopt+1;%复合型进化时取样的数目
nspl=npg;%每个复合型进化的迭代数
mings=ngs;%进化过程复合型的最小数目
npt=npg*ngs;%复合型总的点数
%计算上下界之差
bound=bu-bl;
%初始化种群
rand('seed',iseed);
x=zeros(npt,nopt);
for i=1:npt
    x(i,:)=int32(bl+rand(1,nopt).*bound);
end
%是否包含x0
if iniflg==1
    x(1,:)=x0;
end
%计算函数值
nloop=0;
icall=0;
for i=1:npt
    xf(i)=func(x(i,:));
    icall=icall+1;
end
f0=xf(1);
%按函数值不减排序
[xf,idx]=sort(xf);
x=x(idx,:);
%
%记录最好最差点
bestx=x(1,:);
worstx=x(npt,:);
bestfx=xf(1);
worstfx=xf(npt);
BESTX=bestx;
BESTF=bestfx;
ICALL=icall;
%
%计算标准差
xnstd=std(x);
%计算收敛指标
gnrng=Get_gnrng(x,bound);
%
disp('The Initial Loop: 0');
disp(['BESTF  : ' num2str(bestfx)]);
disp(['BESTX  : ' num2str(bestx)]);
disp(['WORSTF : ' num2str(worstfx)]);
disp(['WORSTX : ' num2str(worstx)]);
disp(' ');
%
if icall>=maxn
    disp('优化搜索因最大数量的限制而终止试验');
    disp(maxn);
    disp('已经超过. 搜索是停在试验数:');
    disp(icall);
    disp('of the initial loop');
end
%
if gnrng<peps
    disp('the population has converged to a prespecified small paremeter space');
end
%
criter=[];
criter_change=1e+5;
while icall<maxn & gnrng>peps & criter_change>pcento
    nloop=nloop+1;
    for igs=1:ngs
        k1=1:npg;
        k2=(k1-1)*ngs+igs;
        cx(k1,:)=x(k2,:);
        cf(k1)=xf(k2);
        
        for loop=1:nspl
            lcs(1)=1;
            for k3=2:nps
                for iter=1:1000
                    lpos=1+floor(npg+0.5-sqrt((npg+0.5)^2-npg*(npg+1)*rand));
                    idx=find(lcs(1:k3-1)==lpos);
                    if isempty(idx)
                        break;
                    end
                end
                lcs(k3)=lpos;
            end
            %对lcs排序
            [lcs,idx1]=sort(lcs);
            %构建单纯性
            s=zeros(nps,nopt);
            s=cx(lcs,:);
            sf=cf(lcs);
            [snew,fnew,icall]=cceua(s,sf,bl,bu,icall,maxn);
            %代替最差的点
            s(nps,:)=snew;
            sf(nps)=fnew;
            %把单纯性放回到复合型
            cx(lcs,:)=s;
            cf(lcs)=sf;
            %对复合型排序
            [cf,idx2]=sort(cf);
            cx=cx(idx2,:);
        end
        %把复合型放回到种群中
        x(k2,:)=cx(k1,:);
        xf(k2)=cf(k1);
    end
   %混合复合型排序
   [xf,idx3]=sort(xf);
   x=x(idx3,:);
   PX=x;
   PF=xf;
   %记录最好最差点
   bestx=x(1,:);
   worstx=x(npt,:);
   bestfx=xf(1);
   worstfx=xf(npt);
   BESTX=[BESTX;bestx];
   BESTF=[BESTF;bestfx];
   ICALL=[ICALL;icall];
   %计算标准差
   xnstd=std(x);
   %计算收敛指标
   gnrng=Get_gnrng(x,bound);
   
   disp(['Evolution Loop: ' num2str(nloop) '  - Trial - ' num2str(icall)]);
   disp(['BESTF  : ' num2str(bestfx)]);
   disp(['BESTX  : ' num2str(bestx)]);
   disp(['WORSTF : ' num2str(worstfx)]);
   disp(['WORSTX : ' num2str(worstx)]);
   disp(' ');
   
   % Check for convergency;
   if icall > maxn
       disp(['已达到最大调用次数：' num2str(maxn)]);
   end
   if gnrng < peps
       disp('种群已进化到了一个预料中的小范围');
   end
   
   criter=[criter;bestfx];
   if (nloop >= kstop)
       criter_change=abs(criter(nloop)-criter(nloop-kstop+1))*100;
       criter_change=criter_change/mean(abs(criter(nloop-kstop+1:nloop)));
       if criter_change < pcento
           disp(['最优点已经非常接近，循环次数为：' num2str(kstop) ' 较阈值：' num2str(pcento) '%']);
       end
   end
   % End of the Outer Loops
end

disp(['计算次数: ' num2str(icall)]);
disp(['归一化几何指数： ' num2str(gnrng)]);
disp(['最好点较上个点进化程度为 ：' num2str(criter_change) '%']);

% END of Subroutine sceua
return;    

function gnrng=Get_gnrng(x,bound)
gnrng=exp(mean(log((max(x)-min(x))./bound)));
return;

function f=func(x)
disp(['调整参数为：' num2str(x)]);
LAI=input('模型反演LAI值：');
RSLAI=[1.42 2.09];
f=mean(abs(LAI-RSLAI));
disp(['函数值为：' f]);
return;

function [snew,fnew,icall]=cceua(s,sf,bl,bu,icall,maxn)
[n,m]=size(s);
alpha=1.0;
beta=0.5;
%记录最好与最差点
bestx=s(1,:);
worstx=s(n,:);
bestfx=sf(1);
worstfx=sf(n);
%计算形心（去掉最差点）
ce=mean(s(1:n-1,:));
%计算反射点
snew=int32(ce+alpha*(ce-worstx));
%检查是否出界
ibound=0;
s1=snew-bl;idx=find(s1<0);if ~isempty(idx) ibound=1;end
s1=snew-bu;idx=find(s1>0);if ~isempty(idx) ibound=2;end
%
if ibound >=1
    snew=int32(bl+rand(1,m).*(bu-bl));
end
fnew=func(snew);
icall=icall+1;
%如果反射点失败，尝试收缩点
if fnew>worstfx
    snew=int32(worstx+beta*(ce-worstx));
    fnew=func(snew);
    icall=icall+1;
    %使用随机点
    if fnew>worstfx
        snew=int32(bl+rand(1,m).*(bu-bl));
        fnew=func(snew);
        icall=icall+1;
    end
end
return;