...
 
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;
warn.order_previous = order_previous';
warn.order_first = order_first';
warn.warning = 'the matchings are not the same';
warningCell{warningCnt} = warn;
warningCnt = warningCnt+1;
% k
% order_previous'
% order_first'
% warning('the matchings are not the same');
end
eigVal_previous = eigVals(order_previous);
eigValAll_previous(k,:) = eigVal_previous;
eigVec_previous = eigVec(:,order_previous);
eigValAll_first(k,:) = eigVal_first;
for v1=1:dim
if si(v1)==-1
eigVec_previous(:,v1)= -1*eigVec_previous(:,v1);
end
end
% order_previous'
% order_first'
end
end
%divide by numLast --> get mean
eigVec_previous = eigVec_m_previous./numLast;
eigVals_previous = eigVal_m_previous./numLast;
% eigVal_m_previous
% eigVal_m_first
eigVec_first = eigVec_m_first./numLast;
eigVals_first = eigVal_m_first./numLast;
%normalize eigvectors v/norm(v)
eigVec_previous = eigVec_previous./norm(eigVec_previous); % for k=1:dim eigVec2(:,k) = eigVec(:,k)./norm(eigVec(:,k));end
eigVec_first = eigVec_first./norm(eigVec_first);
%get covariance
% previous
C_previous = eigVec_previous *diag(eigVals_previous)* inv(eigVec_previous)
%symmetrize C
C_previous_sym = 0.5*(C_previous + C_previous');
% asymmetric part
C_previous_asym = 0.5*(C_previous - C_previous');
% first
C_first = eigVec_first * diag(eigVals_first) * inv(eigVec_first);
%symmetrize C
C_first_sym = 0.5*(C_first + C_first');
% asymmetric part
C_first_asym = 0.5*(C_first - C_first');
%
% r_norm_previous = nthroot(det(C_previous),dim*2);
% r_norm_first = nthroot(det(C_first),dim*2);
%
% Vol_previous = Vol_lp(dim,r_norm_previous,p_test)*out.P_empAll;
% Vol_first = Vol_lp(dim,r_norm_first,p_test)*out.P_empAll;
out_averageCov.Vec_similar_previous = Vec_similar_previous;
out_averageCov.Vec_similar_first = Vec_similar_first;
out_averageCov.Vol_test_all = Vol_test_all;
out_averageCov.C_previous = C_previous_sym;
out_averageCov.C_first = C_first_sym;
out_averageCov.C_previous_asym = C_previous_asym;
out_averageCov.C_first_asym = C_first_asym;
out_averageCov.p_test = p_test;
out_averageCov.dim = dim;
out_averageCov.eigValAll_first = eigValAll_first;
out_averageCov.eigValAll_previous = eigValAll_previous;
% if warningCnt > 1
% out_warningCell = warningCell{1:(warningCnt-1),:};
% end
warningCell = warningCell(1:(warningCnt-1));
out_averageCov.warningCell = warningCell;
\ No newline at end of file
function [h,x,y,z]=error_ellipse(varargin)
% ERROR_ELLIPSE - plot an error ellipse, or ellipsoid, defining confidence region
% ERROR_ELLIPSE(C22) - Given a 2x2 covariance matrix, plot the
% associated error ellipse, at the origin. It returns a graphics handle
% of the ellipse that was drawn.
%
% ERROR_ELLIPSE(C33) - Given a 3x3 covariance matrix, plot the
% associated error ellipsoid, at the origin, as well as its projections
% onto the three axes. Returns a vector of 4 graphics handles, for the
% three ellipses (in the X-Y, Y-Z, and Z-X planes, respectively) and for
% the ellipsoid.
%
% ERROR_ELLIPSE(C,MU) - Plot the ellipse, or ellipsoid, centered at MU,
% a vector whose length should match that of C (which is 2x2 or 3x3).
%