% evolve a Markov chain or discrete dynamical system
% barnett 8/11/06

% MARKOV CHAIN ............................................................
A = [0 1/2 1/2;1 0 1/2; 0 1/2 0];      % abc web
%A = [0 1/2 0; 1 0 0; 0 1/2 1];      % c no outlinks
%A = [0 1/2 1/3; 1 0 1/3; 0 1/2 1/3];      % c no outlinks, fixed via rand jump
%A = [0 0 1/3; 1 0 1/3; 0 1 1/3];      % a no inlinks, c no outlinks

n = size(A,1);
[V,D] = eig(A);                 % diagonalize
diag(D)                         % show eigenvalues of A
v = null(A-eye(n), 'r');        % eigenvector for eigenvalue 1 via nullspace
v / sum(v)                      % predict steady state vector

% evolve it while displaying 3 ways...
T = 20;
x = zeros(n,T+1);
x(1,1) = 1;                     % set initial col vec
figure; set(gcf, 'position', [300 100 400 800]);
for k=1:T
  subplot(4,1,1); plot(1:k, x(:,1:k), '+-'); axis([1 T 0 1]); xlabel('k');
  ylabel('x_j(k)');
  subplot(4,1,2); imagesc(x(:,1:T)); caxis([0 1]); colormap(gray(256));
  xlabel('k'); ylabel('state');
  subplot(4,1,[3 4]); plot(x(1,1:k), x(2,1:k), 'o-'); axis([0 1 0 1]);
  xlabel('x_1(k)'); ylabel('x_2(k)'); set(gca, 'xgrid', 'on', 'ygrid', 'on');
  drawnow; pause(0.2);   % make plotting as it happens, and slow it down
  x(:,k+1) = A*x(:,k);                 % calculate
end
x(:,end)


% DISCRETE LINEAR DYNAMICAL SYSTEM ..........................................

%A = [-5 2; -12 5];   % WS A and E.
%A = [.5 .4; -.104 1.1];             % owls example
%P = [1 1;1 -1]; D = [2 0; 0 1/2]; A = P*D*inv(P);   % worksheet
%P = [1 1;1 -1]; D = [0.9 0; 0 0.8]; A = P*D*inv(P);   % attractor
A = [.5 .2; .5 .8];
eig(A)

n = size(A,1);
% evolve it while displaying 2 ways...
T = 50;
x = zeros(n,T+1);
%x(:,1) = [1;.5];                     % owls initial vector
%x(:,1) = [1;1];                     % WS initial vector
x(:,1) = [.7;.3];                     % MC
figure;
for k=1:T
  subplot(2,1,1); plot(1:k, x(:,1:k), '+-');
  axis([1 T 0 1]); %axis([1 T -1 2]);
  xlabel('k'); ylabel('x_j(k)');
  subplot(2,1,2); plot(x(1,1:k), x(2,1:k), 'o-'); axis equal;%axis([0 2 0 2]); 
  xlabel('x_1(k)'); ylabel('x_2(k)'); set(gca, 'xgrid', 'on', 'ygrid', 'on');
  drawnow; pause(0.05);   % make plotting as it happens, and slow it down
  x(:,k+1) = A*x(:,k);                 % calculate
end
% add the eigenvectors to the second plot... 
[V,D] = eig(A);                     % diagonalize
hold on; plot([V(1,:); -V(1,:)], [V(2,:); -V(2,:)], 'r-*');

% end
