function DS = ON_2DS(O, iterations)

n = size(O,1);
if(~exist('iterations','var'))
    iterations = 1000;
end

DS = zeros(size(O));
for i=1:iterations
    x = randn(n,1);
    %x=x./norm(x);
    y = O*x;

    [~,indx] = sort(x,'ascend');
    [~,indy] = sort(y,'ascend');

    ind = sub2ind(size(O), indy, indx);
    P = zeros(size(O));
    P(ind) = 1;
    DS = DS + P;
end

DS = DS./iterations;

end