...
 
Commits (2)

Too many changes to show.

To preserve performance only 1000 of 1000+ files are displayed.

File added
function xmin=EAforStartPoint % (mu/mu_w, lambda)-CMA-ES
% -------------------- Initialization --------------------------------
% User defined input parameters (need to be edited)
strfitnessfct = 'evaluateRegions'; % name of objective/fitness function
N = 3; % number of objective variables/problem dimension
xmean = rand(N,1); % objective variables initial point
sigma = 0.3; % coordinate wise standard deviation (step size)
stopfitness = 0; % stop if fitness < stopfitness (minimization)
stopeval = 1e3*N^2; % stop after stopeval number of function evaluations
% Strategy parameter setting: Selection
lambda = 4+floor(3*log(N)); % population size, offspring number
mu = lambda/2; % number of parents/points for recombination
weights = log(mu+1/2)-log(1:mu)'; % muXone array for weighted recombination
mu = floor(mu);
weights = weights/sum(weights); % normalize recombination weights array
mueff=sum(weights)^2/sum(weights.^2); % variance-effectiveness of sum w_i x_i
% Strategy parameter setting: Adaptation
cc = (4+mueff/N) / (N+4 + 2*mueff/N); % time constant for cumulation for C
cs = (mueff+2) / (N+mueff+5); % t-const for cumulation for sigma control
c1 = 2 / ((N+1.3)^2+mueff); % learning rate for rank-one update of C
cmu = min(1-c1, 2 * (mueff-2+1/mueff) / ((N+2)^2+mueff)); % and for rank-mu update
damps = 1 + 2*max(0, sqrt((mueff-1)/(N+1))-1) + cs; % damping for sigma
% usually close to 1
% Initialize dynamic (internal) strategy parameters and constants
pc = zeros(N,1); ps = zeros(N,1); % evolution paths for C and sigma
B = eye(N,N); % B defines the coordinate system
D = ones(N,1); % diagonal D defines the scaling
C = B * diag(D.^2) * B'; % covariance matrix C
invsqrtC = B * diag(D.^-1) * B'; % C^-1/2
eigeneval = 0; % track update of B and D
chiN=N^0.5*(1-1/(4*N)+1/(21*N^2)); % expectation of
% ||N(0,I)|| == norm(randn(N,1))
% -------------------- Generation Loop --------------------------------
counteval = 0; % the next 40 lines contain the 20 lines of interesting code
while counteval < stopeval
disp('pitsi');
disp(counteval);
% Generate and evaluate lambda offspring
for k=1:lambda
arx(:,k) = xmean + sigma * B * (D .* randn(N,1)); % m + sig * Normal(0,C) for i=1:N
if(arx(1,k)<0)
arx(1,k)=0;
end
if(arx(1,k)>0.3)
arx(1,k)=0.3;
end
if(arx(2,k)<0)
arx(2,k)=0;
end
if(arx(2,k)>1)
arx(2,k)=1;
end
if(arx(3,k)<0)
arx(3,k)=0;
end
arfitness(k) = feval(strfitnessfct, arx(:,k)); % objective function call
counteval = counteval+1;
end
% Sort by fitness and compute weighted mean into xmean
[arfitness, arindex] = sort(arfitness); % minimization
xold = xmean;
xmean = arx(:,arindex(1:mu))*weights; % recombination, new mean value
% Cumulation: Update evolution paths
ps = (1-cs)*ps ...
+ sqrt(cs*(2-cs)*mueff) * invsqrtC * (xmean-xold) / sigma;
hsig = norm(ps)/sqrt(1-(1-cs)^(2*counteval/lambda))/chiN < 1.4 + 2/(N+1);
pc = (1-cc)*pc ...
+ hsig * sqrt(cc*(2-cc)*mueff) * (xmean-xold) / sigma;
% Adapt covariance matrix C
artmp = (1/sigma) * (arx(:,arindex(1:mu))-repmat(xold,1,mu));
C = (1-c1-cmu) * C ... % regard old matrix
+ c1 * (pc*pc' ... % plus rank one update
+ (1-hsig) * cc*(2-cc) * C) ... % minor correction if hsig==0
+ cmu * artmp * diag(weights) * artmp'; % plus rank mu update
% Adapt step size sigma
sigma = sigma * exp((cs/damps)*(norm(ps)/chiN - 1));
% Decomposition of C into B*diag(D.^2)*B' (diagonalization)
if counteval - eigeneval > lambda/(c1+cmu)/N/10 % to achieve O(N^2)
eigeneval = counteval;
C = triu(C) + triu(C,1)'; % enforce symmetry
[B,D] = eig(C); % eigen decomposition, B==normalized eigenvectors
D = sqrt(diag(D)); % D is a vector of standard deviations now
invsqrtC = B * diag(D.^-1) * B';
end
% Break, if fitness is good enough or condition exceeds 1e14, better termination methods are advisable
if arfitness(1) <= stopfitness || max(D) > 1e7 * min(D)
break;
end
end % while, end generation loop
xmin = arx(:, arindex(1)); % Return best point of last iteration.
% Notice that xmean is expected to be even
% better.
% ---------------------------------------------------------------
function f=frosenbrock(x)
if size(x,1) < 2 ,error('dimension must be greater one'); end
f = 100*sum((x(1:end-1).^2 - x(2:end)).^2) + sum((x(1:end-1)-1).^2);
\ No newline at end of file
This diff is collapsed.
function [Vol] = Vol_lp(N,r,p)
% N - dimension
% r - radius
% p - p-norm
if p > 100
Vol = (2*r).^N;
else
Vol = (2*gamma(1/p+1).*r).^N/(gamma(N/p+1));
end
\ No newline at end of file
<
function [out_averageCov,warningCell]= averageCov(out,numLast)
p_test = out.opts.pn;
dim = size(out.xAcc,2);
% save estimated volume from MC at each iteration (of the last numLast iter)
Vol_test_all = nan(numLast,1);
% save how close eigenvectors are to first set of eigenvectors
Vec_similar_first = nan(numLast-1,dim);
Vec_similar_previous = nan(numLast-1,dim);
% save average eigenvalues and eigenvectors
eigVal_m_previous = nan(dim,1);
eigVec_m_previous = nan(dim);
eigVal_m_first = nan(dim,1);
eigVec_m_first = nan(dim);
eigValAll_first = nan(numLast,dim);
eigValAll_previous = nan(numLast,dim);
warningCell=cell(numLast,1);
warningCnt = 1;
for k=1:numLast
%numLast
%size(out.QCell)
%get MC Covariance in each step
Q = out.QCell{end-numLast+k};
r = out.rVec(end-numLast+k);
C = r^2*(Q*Q');
%estimated Volume from MC
Vol_test_all(k) = Vol_lp(dim,r,p_test)*out.P_empVecAll(end-numLast+k);
%eigen decomposition
[eigVec,eigVals] = eig(C);
if k==1
%to save average eigenvalues and eigenvectors
eigVal_m_previous=diag(eigVals);
eigVec_m_previous=eigVec;
eigVal_m_first = diag(eigVals);
eigVec_m_first = eigVec;
%eigenvalues and eigenvectors from previous iteration
eigVal_previous = diag(eigVals);
eigVec_previous = eigVec;
%eigenvalues and eigenvectors from first iteration
eigVal_first = diag(eigVals);
eigVec_first = eigVec;
eigValAll_previous(k,:) = eigVal_previous;
eigValAll_first(k,:) = eigVal_first;
else
check_allUsed_previous = zeros(dim,1);
check_allUsed_first = zeros(dim,1);
order_previous = zeros(dim,1);
order_first = zeros(dim,1);
%signs
si = ones(dim,1);
for v1=1:dim
%for each eigenvector, find corresponding eigenvector from
%previous iteration and from .
eigVec1 = eigVec(:,v1);
%get vectors that are most alike from previous iteration
[maxi_previous,idx_previous]=max(abs(eigVec1'*eigVec_previous));
check_allUsed_previous(idx_previous) = check_allUsed_previous(idx_previous)+1;
order_previous(v1)=idx_previous;
%get vectors that are most alike from first iteration
[maxi_first,idx_first] = max(abs(eigVec1'*eigVec_first));
check_allUsed_first(idx_first) = check_allUsed_first(idx_first)+1;
order_first(v1) = idx_first;
%add eigenvalues
diag_eigVals = diag(eigVals);
eigVal_m_first(idx_first) = eigVal_m_first(idx_first)+diag_eigVals(v1);
eigVal_m_previous(idx_previous) = eigVal_m_previous(idx_previous)+diag_eigVals(v1);
% % check that vectors are pointing in same direction
% % check angle between them
% % multiply eigVec1 * -1 if angle too big
% angle_v = acos(eigVec1'*eigVec_old(:,idx));
% if angle_v >= pi/4
% eigVec1 = eigVec1*(-1);
% disp('here');
% end
% check that vectors are pointing in same direction
projection_previous = eigVec1'*eigVec_previous(:,idx_previous);
if projection_previous < 0
eigVec1_previous = eigVec1*(-1);
si(v1)=-1;
else
eigVec1_previous = eigVec1;
end
projection_first = eigVec1'*eigVec_first(:,idx_first);
if projection_first < 0
eigVec1_first = eigVec1*(-1);
else
eigVec1_first = eigVec1;
end
% add eigenvectors
eigVec_m_previous(:,idx_previous)=eigVec_m_previous(:,idx_previous)+eigVec1_previous;
eigVec_m_first(:,idx_first)=eigVec_m_first(:,idx_first)+eigVec1_first;
%vec similar
Vec_similar_previous(k-1,idx_previous) = maxi_previous;
Vec_similar_first(k-1,idx_first) = maxi_first;
end
if ~isempty(find(check_allUsed_previous==0))% any(iszero(check_allUsed_previous))
k
check_allUsed_previous
error('not all eigenvectors matched to');
end
if ~isempty(find(check_allUsed_first==0))%any(iszero(check_allUsed_first))
k
warn.k=k;
warn.check_allUsed_first=check_allUsed_first;
warn.order_previous = order_previous;
warn.order_first=order_first;
warn.warning='when matching to the first iteration, at least one eigenvector is not matched to';
% k
% check_allUsed_first
% order_previous
% order_first
% warning('when matching to the first iteration, at least one eigenvector is not matched to');
warningCell{warningCnt} = warn;
warningCnt = warningCnt+1;
end
if ~all(order_previous == order_first)
k
warn.k=k;