function [s ranks rankDelta] = noisemeasurement(arf1, arf2, lamreev, theta, cutlimit) % function [s ranks rankDelta] = noisemeasurement(arf1, arf2, lamreev, theta) % % Input: % arf1, arf2 : two arrays of function values. arf1 is of size 1xlambda, % arf2 may be of size 1xlamreev or 1xlambda. The first lamreev values % in arf2 are (re-)evaluations of the respective solutions, i.e. % arf1(1) and arf2(1) are two evaluations of "the first" solution. % lamreev: number of reevaluated individuals in arf2 % theta : parameter theta for the rank change limit, between 0 and 1, % typically between 0.2 and 0.7. % cutlimit (optional): output s is computed as a mean of rankchange minus % threshold, where rankchange is <=2*(lambda-1). cutlimit limits % abs(rankchange minus threshold) in this calculation to cutlimit. % cutlimit=1 evaluates basically the sign only. cutlimit=2 could be % the rank change with one solution (both evaluations of it). % % Output: % s : noise measurement, s>0 means the noise measure is above the % acceptance threshold % ranks : 2xlambda array, corresponding to [arf1; arf2], of ranks % of arf1 and arf2 in the set [arf1 arf2], values are in [1:2*lambda] % rankDelta: 1xlambda array of rank movements of arf2 compared to % arf1. rankDelta(i) agrees with the number of values from % the set [arf1 arf2] that lie between arf1(i) and arf2(i). % % Note: equal function values might lead to somewhat spurious results. % For this case a revision is advisable. %%% verify input argument sizes if size(arf1,1) ~= 1 error('arf1 must be an 1xlambda array'); elseif size(arf2,1) ~= 1 error('arf2 must be an 1xsomething array'); elseif size(arf1,2) < size(arf2,2) % not really necessary, but saver error('arf2 must not be smaller than arf1 in length'); end lam = size(arf1, 2); if size(arf1,2) ~= size(arf2,2) arf2(end+1:lam) = arf1((size(arf2,2)+1):lam); end if nargin < 5 cutlimit = inf; end %%% capture unusual values if any(diff(arf1) == 0) % this will presumably interpreted as rank change, because % sort(ones(...)) returns 1,2,3,... warning([num2str(sum(diff(arf1)==0)) ' equal function values']); end %%% compute rank changes into rankDelta % compute ranks [ignore, idx] = sort([arf1 arf2]); [ignore, ranks] = sort(idx); ranks = reshape(ranks, lam, 2)'; rankDelta = ranks(1,:) - ranks(2,:) - sign(ranks(1,:) - ranks(2,:)); %%% compute rank change limits using both ranks(1,...) and ranks(2,...) for i = 1:lamreev sumlim(i) = ... 0.5 * (... myprctile(abs((1:2*lam-1) - (ranks(1,i) - (ranks(1,i)>ranks(2,i)))), ... theta*50) ... + myprctile(abs((1:2*lam-1) - (ranks(2,i) - (ranks(2,i)>ranks(1,i)))), ... theta*50)); end %%% compute measurement %s = abs(rankDelta(1:lamreev)) - sumlim; % lives roughly in 0..2*lambda % max: 1 rankchange in 2*lambda is always fine s = abs(rankDelta(1:lamreev)) - max(1, sumlim); % lives roughly in 0..2*lambda % cut-off limit idx = abs(s) > cutlimit; s(idx) = sign(s(idx)) * cutlimit; s = mean(s); %%% improvements to the original publication: % 1) (DONE) s = mean(sign(2*abs(rankDelta(1:lamreev)) - sumlim)); (using % the sign is better in case of outliers), and take a mean in the % generation sequence, e.g. such that lamreev\approx10. Nasty test: % say 40% outliers inside a ball around the optimum, such that the % sigma-increase is not sufficient to escape (for lamreev=2), but % that the selection is enough disturbed to get stuck. % % 2) (DONE) To avoid that 1/2 exchange increases sigma: sumlim(i) = max(1, % prctile(...)) + max(1,prctile(...)); the effect is quite % remarkable (much close to the optimum, but with cum=1?) but it % is not clear whether it is actually preferable. % % see also testrank function res = myprctile(inar, perc, idx) % % Computes the percentiles in vector perc from vector inar % returns vector with length(res)==length(perc) % idx: optional index-array indicating sorted order % N = length(inar); flgtranspose = 0; % sizes if size(perc,1) > 1 perc = perc'; flgtranspose = 1; if size(perc,1) > 1 error('perc must not be a matrix'); end end if size(inar, 1) > 1 & size(inar,2) > 1 error('data inar must not be a matrix'); end % sort inar if nargin < 3 || isempty(idx) [sar idx] = sort(inar); else sar = inar(idx); end res = []; for p = perc if p <= 100*(0.5/N) res(end+1) = sar(1); elseif p >= 100*((N-0.5)/N) res(end+1) = sar(N); else % find largest index smaller than required percentile availablepercentiles = 100*((1:N)-0.5)/N; i = max(find(p > availablepercentiles)); % interpolate linearly res(end+1) = sar(i) ... + (sar(i+1)-sar(i))*(p - availablepercentiles(i)) ... / (availablepercentiles(i+1) - availablepercentiles(i)); end end if flgtranspose res = res'; end % Hansen, N., A.S.P. Niederberger, L. Guzzella and P. Koumoutsakos % (2009?). A Method for Handling Uncertainty in Evolutionary % Optimization with an Application to Feedback Control of % Combustion. To appear in IEEE Transactions on Evolutionary % Computation.