% The following MATLAB .m file finds the value of the second highest
% eigenvalue in absolute value of the matrix P from the 'frog example'.
% A lower second highest eigenvalue generally means the matrix converges
% more quickly to stationarity.

m = 9;       % P(x,y) = 1/(2m+1) for all x, y within m of each other
nodes = 200; % P is a nodes*nodes matrix
beta = 10;   % target distribution is proportional to (dist(x,1)+1)^-b
             % where dist(x,1) is the number of nodes to get from x to 1

% the following two loops calculate the stationary (target) distribution
% of the matrix P
denom = 0;
for i = 1:nodes
  min = i-1;
  if nodes-i+1 < min
    min = nodes-i+1;  % distance from the current node to 1
  end
  denom = denom + (min + 1)^(-beta);
end
for i = 1:nodes
  min = i-1;
  if nodes-i+1 < min
    min = nodes-i+1;
  end
  pi(1, i) = (min + 1)^(-beta)/denom;  % stationary distribution of P
end

for i = 1:nodes
   P(i, i) = 1/(2*m+1);
   for j = 1:nodes
      if pi(1, j) >= pi(1, i)
         alpha = 1;  % alpha is the probability that we jump from the
                     % current node, i, to the proposed node, j
      else
         alpha = pi(1, j)/pi(1, i);
      end

      if 0 < abs(i-j) & abs(i-j) < (m+1)
         P(i, j) = alpha/(2*m+1);
         P(i, i) = P(i, i) + (1-alpha)/(2*m+1);
      elseif abs(i-j) >= (nodes-m)  % the highest nodes and the lowest
                                    % nodes have edges between them too
         P(i, j) = alpha/(2*m+1);
         P(i, i) = P(i, i) + (1-alpha)/(2*m+1);
      elseif (m+1) <= abs(i-j) & abs(i-j) < (nodes-m)
         P(i, j) = 0;
      end
   end
end

y = eig(P);  % finds a vector of the eigenvalues; the following loop
             % isolates the second highest one in absolute value
max = 0;
max2 = 0;
for k = 1:nodes
   if abs(y(k, 1)) > max
      max2 = max;
      max = abs(y(k, 1));
   elseif abs(y(k, 1))>max2 & max>= abs(y(k, 1))
      max2 = abs(y(k, 1));
   end
end
max2


% the following can be added on to the above program to find bounds
% on variation distance using all the eigenvalues; a slight modification
% can be made to see bounds on L2(pi) distance
iters = 40;
Q = P'; % we want left eigenvectors of P which will be right eigenvectors
        % of Q
sum = 0;
[evecs, evals] = eig(Q);
start = zeros(nodes, 1);
start(1,1) = 1;
am = evecs\start;  % finding a sub m for all m
for i = 1:nodes
  for j = 2:nodes
    sum = sum + abs(am(j,1)*evecs(i,j)) * abs(evals(j,j))^iters;
  end
end
sum = sum/2;
sum  % bound on variation distance