  function P4Grade()
% Checks out pLimits
% 10 points
% Roughly half credit for methods that use svd or eig instead of eigs

clc
% Use n=5 to see what is going on
Q = randQ(5,1);
F = expm(Q*20);
% We get something like this
%        0.0998    0.0998    0.0998    0.0998    0.0998
%        0.1656    0.1656    0.1656    0.1656    0.1656
%  F =   0.2278    0.2278    0.2278    0.2278    0.2278
%        0.2472    0.2472    0.2472    0.2472    0.2472
%        0.2596    0.2596    0.2596    0.2596    0.2596
%
% Note that F = v*ones(n,1)' and that diag(F) is what we want.
% What is v?
% 
% To see what is going on, suppose Q is diagonalizable:
%     Q = Xdiag(d)inv(X) = d(1)X(:,1)Y(:,1)' + ...+ d(n)X(:,n)Y(:,n)'
% where inv(X) = Y'. It follows that
%     exp(Qt) = Xdiag(d)inv(X) = exp(d(1)t)X(:,1)Y(:,1)' + ...+ exp(d(n)t)X(:,n)Y(:,n)'
% 
% Since this converges to a rank 1, all the d(i) except one converges to
% zero. We know that Q has one eigenvalu equal to zero, so that is the 
% term that "sticks around". Thus, if d(1) = 0 then
%       exp(Qt) \approx X(:,1)Y(:,1)'.
% QX(:,1) = 0 and Y(:,1)'Q = 0. Because of the zero column sumns, Y(:,1) =
% ones(n,1). 
% Bottom line: compute the (right) null vector for Q. 
% We will want to avoid an O(n^3) factorization since we just need that
% one vector. Let's apply power method to a shifted version of Q (alfa*I+Q)
% If alfa is positive enough, then it will be the (unique) dominant eigenvalue
% and power methoid is OK. alfa = norm(Q,1) works. And instead of power method
% lets call eigs.

% A small example for correctness

n = 10;
Q = randQ(n,10);
pYOU = pLimits(Q);
err = 1; t = 0;
fprintf('Convergence Example, n = %1d\n',n)
fprintf('    t   || pYOU - diag(expm(Q*t)) ||\n--------------------------------\n')
while err>=10^-12 && t <=100
    err = norm(diag(expm(t*Q))-pYOU,1);
    fprintf('   %2.0f          %10.3e\n',t,err)
    t = t+1;
end

 
% Large example for efficiency
n = 1000;
Q = randQ(n,1);
tic
pYOU = pLimits(Q);
tYOU = toc;
tic
pCVL = pLimitsCVL(Q);
tCVL = toc;
fprintf('\nn = %4d  tYOU/tCVL = %6.3f\n',n,tYOU/tCVL)

  function p = pLimitsCVL(Q)
% Q is an intensity matrix and p = lim diag(expm(Qt))
[n,n] = size(Q);
alfa = norm(Q,1);
A = alfa*eye(n,n)+Q;
[p,e] = eigs(A,1);
% Normalize...
p = p/sum(p);


  function Q = randQ(n,tau)
% Random intensity matrix.
% n is apositive integer and tau is a positive float. The bigger the value
% of tau, the more variation along Q's diagonal: diag(Q) = -tau*rand(n,1),
d = tau*rand(n,1);
Q = zeros(n,n);
for j=1:n
    v = rand(n-1,1); 
    v = d(j)*v/sum(v);
    Q([1:j-1 j+1:n],j) = d(j)*v/sum(v);
    Q(j,j) = -d(j);
end
    
