function v = invariantSIM(P,N)
%
% v = invariantSIM(P,N)
%
n = length(P);
i = randdisc(ones(1,n)/n);
v = zeros(1,n);
for k = 1:N
  i = randdisc(P(i,:));
  v(i) = v(i)+1;
end
v = v/N;
