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.