amoeba.asv

标准单纯性算法程序

function bestfit=amoeba(funk,xguess,iterations)

p=size(xguess,2);
beta=1;				 % Reflection coefficient (standard choice)
gamma=1;		     % Expansion coefficient (standard choice)
alpha=0.5;			 % Contraction coefficient (standard choice)
[M,N] = size(xguess); %N will be the number of dimensions.
if M~=1 
   error('xguess is not a row vector!')
end

%Define the N+1 'feet' of the amoeba
basis = eye(N);
feet = [xguess;ones(N,1)*xguess+basis]';
P(:,:,1) = feet;
for j=1:iterations
    j

	% First, compute the values of all the vertices
	
	for i=1:p+1
		
J(i,j)=feval(funk,[P(1,i,j);P(2,i,j)]);

   end

% Find the worst and best vertex and their indices

[Jmax(j),maxvertex(j)]=max(J(:,j));
[Jmin(j),minvertex(j)]=min(J(:,j));

thetamax(:,j)=P(:,maxvertex(j),j);
thetamin(:,j)=P(:,minvertex(j),j);

% Find the centroid of all the points, except the worst one

thetacent(:,j)=0*ones(p,1);  % Initialize it

for i=1:p+1
	
  thetacent(:,j)=thetacent(:,j)+P(:,i,j);

end

thetacent(:,j)=(1/p)*(thetacent(:,j)-thetamax(:,j));                       %求质心

% Reflection point computation step:

thetaref(:,j)=thetacent(:,j)+beta*(thetacent(:,j)-thetamax(:,j));          %反射

Jref(j)=feval(funk,[thetaref(1,j);thetaref(2,j)]);

% Next compute some values for the inequality tests in the reflection step

temp1=J(maxvertex(j),j); % Save the value of the cost at the worst vertex
J(maxvertex(j),j)=-Inf;  % Replaces the max element with a small element 
										% (for this problem then this element will not be 
										% chosen as the max)
[Jmaxwm(j),maxvertexwm(j)]=max(J(:,j)); % Finds the second largest element                 %求次大值
J(maxvertex(j),j)=temp1; % Returns the matrix to how it was

% Next perform the inequality tests in the reflection step and then each 
% of the corresponding steps that they indicate.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For the first test in the reflection step

if Jmin(j)>Jref(j),
	% Then perform expansion
	
	thetaexp(:,j)=thetaref(:,j)+gamma*(thetaref(:,j)-thetacent(:,j));

	% Next, we have to compute Jexp, the cost at the expansion point
	
	Jexp(j)=feval(funk,[thetaexp(1,j);thetaexp(2,j)]);
 
   % Then, we attempt expansion
   
	if Jexp(j)<Jref(j),
		thetanew(:,j)=thetaexp(:,j);
		Jnew(j)=Jexp(j);
	else
		thetanew(:,j)=thetaref(:,j);
		Jnew(j)=Jref(j);
	end

	
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For the second test in reflection step

if Jmaxwm(j)>Jref(j) & Jref(j)>=Jmin(j),   
	% Reflection point has an intermediate value
	% so use the "use reflection step"
	
	thetanew(:,j)=thetaref(:,j);
	Jnew(j)=Jref(j);
	
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For the third test in reflection step

if Jref(j)>=Jmaxwm(j), 
	% The reflection point is not good
	
	% We perform contraction step
	
	if Jmax(j)<=Jref(j)
		thetanew(:,j)=alpha*thetamax(:,j)+(1-alpha)*thetacent(:,j);
	else
		thetanew(:,j)=alpha*thetaref(:,j)+(1-alpha)*thetacent(:,j);
	end
	
% We need Jnew for shrinking
	
	Jnew(j)=feval(funk,[thetanew(1,j);thetanew(2,j)]);

end

% Form the new simplex (use shrinking if Jnew in contraction shows 
% that thetanew is worse than the worst vertex):

if Jref(j)>=Jmaxwm(j) & Jnew(j)>Jmax(j), 
				% Tests that came from contraction, and
	            % that the new vertex is not good
		for i=1:p+1
			P(:,i,j+1)=0.5*(P(:,i,j)+thetamin(:,j));  % Shrink towards best vertex
		end
else 
		% Form the simplex in the usual way if you got here via
	    % any step except the the case where shrinking was already used	
	    P(:,:,j+1)=P(:,:,j);
	    P(:,maxvertex(j),j+1)=thetanew(:,j);
end
end % End main loop...
bestfit = P(:,:,end);